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ABSTRACT 

A survey of numerical methods for time dependent partial 
differential equations is presented. The emphasis is on 
practical applications to large scale problems. A discussion 
of new developments in high order methods and moving grids is 
given. The importance of boundary conditions is stressed for 
both internal and external flows. A description of implicit 
methods is presented including generalizations to multidimensions. 
Shocks, aerodynamics, meteorology, plasma physics and combustion 
applications are also briefly described. 
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Introduction 


In the construction of numerical solutions to large 
scale problems one wishes to have a code with many different 
characteristics. Exan^les of these properties include 
efficiency both in terms of computer time and computer 
storage and also ease of use. Frequently a code is con- 
structed to be used by other people who do not have a 
detailed knowledge of the algorithm. Hence, one does not 
want a program which requires the specification of many 
nonphysical parameters or one that requires intervention 
on the part of the user. Lagrangian codes are frequently 
complex and may require operator intervention for rezoning. 
Hence, we shall concentrate on Eulerian methods. 

In contrast to ordinary differential equations it does 
not seem advisable to introduce general packages for time 
dependent partial differential equations. In one space 
dimensions a number of packages exist usually based on an 
ODE solver. These packages are useful when one wants a 
quick answer to a simple problem. However, the programs 
are far from optimal both in terms of computer time and 
computer storage [115]. Hence, for realistic physical 
models with many complicated equations which are to be 
solved many times it is necessary to develop a code for 
each problem. This is especially true for multidimensional 
problems. The range of solutions including both smooth 
and discontinuous flows demands that the algorithm be 
carefully matched with the physics. Depending on many 
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factors one may need either an explicit or implicit scheme. 
Accuracy and geometrical considerations will determine 
whether high or low accuracy methods (in both space and 
time) are more appropriate. 

In the following sections we will discuss in more 
detail many of the factors that influence the choice of a 
scheme. Special attention will be paid to the treatment 
of boundaries. After a general discussion of standard 
boundary treatments attention will be focused on moving 
boundaries. Artificial boundaries to simulate an infinite 
domain create other difficulties which will be analyzed. 

In all these cases the interplay between the physics and 
the numerics will be stressed. Implicit methods and appli- 
cations to specific problems are discussed with special 
emphasis on shocks and on steady state solutions. 


2 . Boundary Treatment 


Since much of this study deals with nonstandard bound- 
ary problems we shall first review the basic theory. For a 
numerical method to be useful we require that it be stable 
and also converge. Hence, small perturbations in the 
problem should give rise to small perturbations in the solu- 
tion. We first consider the model equation 

(2.1) u^ + Au^ = 0 

with A a constant n^n matrix and 0 £ x £ 1. Let A be 
symmetric and have k positive eigenvalues, n-k negative 
eigenvalues. It is then straightforward to show that (2.1) 
is well posed only if k linearly independent conditions are 
imposed at x = 0 and n-k conditions at x = 1. This number 
of conditions is also sufficient for well posedness as long 
as no variables corresponding to characteristic variables 
coming into the boundary, from inside the region, are specified. 

For a numerical method it is necessary to specify the 
correct number of boundary conditions as given by the 
differential equation. In general one requires some method 
to numerically determine the boundary values for the other, 
nonspecified variables. If the boundary treatment is not 
done correctly then errors are generated at the boundary 
which propagate into the domain of integration and create 
instabilities. Especially for nonlinear problems these 
instabilities frequently do not manifest themselves at the 
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boundaries but rather along sonic lines or stagnation 
points. It is a nontrivial task to trace the source of 
- the difficulty to a mistreatment of the boundary [137], 

For simplicity we shall assume that the numerical 
scheme uses information only at the point of interest and 
at its immediate two neighbors at various time levels. Then 
Kreiss has shown [81J, [110] that the stability for a 
scalar equation can be analyzed by assuming a solution of 
the form 

u" = U z^ 

3 

[ and only considering the semi-infinite domain 0 ^ x < 

s 

[ If there are solutions to the interior and boundary differ- 

I 

^ ence schemes with |<| <1 and |z| >1 then the initial- 

i 

[ boundary scheme is unstable. If there are no nontrivial 

solutions with | <' | < 1 and \z\ >1 then the scheme is 
stable. If there are solutions with | ic | =1, \z\ =1 
more care is required. For additional details see the survey 
by Morton [140] . 

For systems of equations or schemes that require more 
than three mesh points at a time level the theory is more 
complicated. For systems, the algebra of solving the equa- 
tions for complex < and z is large and needs to be done 
for each system of equations. Gottlieb, Gunzburger and 
Turkel [72] have shown how the scalar results can be extended 
to systems of equations if one takes into account the 
characteristic variables (see also [28]). If this is not 
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done examples are given to show that troubles can occur. 
Several examples are presented in later chapters. Engquist 
and Smedsaas [55] have a method of lines code that automa- 
tically accounts for the characteristic variables. 

For several of the standard schemes, boundary condi- 
tions have been analyzed by several authors. We shall 
present a brief outline of the results. 

1. Leapfrog method: 

a. Extrapolation in space is unstable. 

b. Extrapolation in space-time is stable. 

c. A one-sided Euler method is stable. 

We again emphasize that these results are true only for a 
scalar equation. As mentioned above the scalar results can 
be extended to systems if the numerical boundary treatment 
is done on the correct characteristic variables rather than 
on the natural variables. The given boundary conditions 
can consist of any combination of variables that yields a 
well posed differential problem. 

2. Lax-Wendroff method: 

This is usually implemented by a two step Richtmyer 
or a two step MacCormack method. The MacCormack method 
is easier to implement as no half points are required and 
boundary conditions can be imposed after the first step. 

The treatment of parabolic terms is also easier using 

the MacCormack scheme. 

a. Space extrapolation at the conclusion of the two 
steps is stable. 

b. A one-sided Euler method at the conclusion of the 
two steps is stable. 


c. Extrapolation at the intermediate step of the 
Richtmyer method is unstable. 

d. Extrapolation or one-sided differences at the 
intermediate step or final step of the MacCormack 
method is stable. 

These boundary treatments are discussed in detail in [71]* 
Method d is particularly easy to implement and will be 
discussed in more detail in the next section where it is 
extended to higher order methods. In many cases one sided 
differences are equivalent to extrapolation of fluxes to 
artificial points exterior tr the region. The choice 
between these options is b' sed on programming ease. 

3. Implicit methods: 

a. Space extrapolation is stable. 

b. The box scheme is stable and very accurate. 

c. Explicit one sided differences can lead to 
stability limits on At ([172]). 

Many numerical tests have confirmed the analysis of 
Kreiss for both simple test cases and complicated problems. 
At present one difficulty is to extend these results, in a 
useful manner, to multidimensional problems. Abarbanel and 
Gottlieb [3] have considered the leapfrog scheme while 
Bayliss (private communication) has analyzed methods based 
on splitting techniques. In practice most of the one 
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dimensional results generalize to the multidimensional 
situation. In the coming chapters we shall consider 
practical generalizations to higher order methods, moving 
boundaries and radiation boundaries. 

The numerical treatment of boundaries for parabolic 
equations is usually simpler than that described above 
since all the variables are prescribed. Some difficulties 
may arise when derivative boundary conditions are given, 
this is especially true for nonrectangular regions [195]. 

More serious difficulties occur when the highest space 
derivatives are multiplied by a small (but fixed) parameter, 
e.g. high Reynolds number flow. Most schemes contain some 
numerical viscosity and so one must ensure that the numeri- 
cal viscosity does not overwhelm the physical viscosity. 

In general this is only true in boundary layers where the 
large gradients enhance the physical viscosity. Therefore, 
when the computations do not resolve the boundary layer it 
is not reasonable to impose parabolic-type boundary condi- 
tions. When one uses a coarse mesh near the boundary one 
is effectively ignoring the viscous effects and only con- 
sidering the inviscid equations. So, with coarse grids near 
the boundary one should use hyperbolic-like boundary condi- 
tions. Failure to do this leads to cell Reynolds number 
restrictions (158). 

To clarify this point we consider a model problem; the 
steady state linearized Burgers equation 

(2.2a) u = Ru , R>0, 
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with boundary conditions 

(2.2b) u(0) = 1 u(l) = 0 . 


The solution is given by 
(2.3) u(x) = (e^ - e^^)/(e^-l) . 


For R large the solution is approximately equal to 1 
everywhere except for a boundary layer near x = 1 . 
Solving (2.2) by central differences we have 


(2. 4a) ''j.i - + '’j-i 


RjIVj+i-Vj.i) with 


(2.4b) R^ = RLx . Vq = 1 , Vj^ = 0 

This has the solution 


(2.5) 


V . 

3 




W = 


2-R^ 

2+R~ 


R, is generally called the cell Reynolds number or the 
Peclet number. For R, larger than 2 , Q is negative and 
so Vj is oscillatory. As R^ increases , Q approaches 
-1 and Vj acquires a large 2Ax oscillation which has 
nothing in common with the analytic solution (see also [158]). 

On« way to avoid this situation is to construct schemes, 
at least near boundaries, which do not have any coll Reynolds 
number restriction [6], [36], [41]. Alternatively one can 

match the interior scheme with an asymptotic expansion for 
the boundary layer [82] . Where feasible stretched grids 
should be used to resolve the boundary layers. When the 


position of the boundary layer is not known in advance or 
is nonuniform along the boundary this is a difficult process. 

The important point is that this oscillation is entirely 
due to faulty boundary treatment, i.e., nonresolution of the 
boundary layers. ?or example, specifying “ 0 at x ■ 1 
rather than u ■ 0 eliminates these oscillations. This is 
a practical method for outflow boundary conditions [84]. 
Alternatively, one can specify combinations of u and higher 
order derivatives at the boundary. As Ax goes to zero this 
combination should reduce to u(l) <*0 and so the scheme con- 
verges for fixed R. As increases the boundary condition 
should approximate some extrapolation formula which is stable 
for the hyperbolic difference approximation. For example, we 
can replace (2.4b) by 

(2.40 ''n ^ ■ “ 


The solution to (2.4) is tlien given by 
(2.6) = AQ^ + 1 - A ; A - ( 1 - (1-Q) ) ‘^ 


By inspection v^ converges to u(x^) as .^x goes to 
zero. As R.^ increases beyond 2, Q becomes negative 
and oscillations appear. However, the factor A in front 
of the oscillatory part becomes small and so the oscillations 
do not disturb the solution. In more practical situations 
the first difference v„ - v„ , can be replaced by higher 
differences for greater accuracy. Similarly, more compli- 
cated weights than simply R can be constructed. For 


10 


multidimensional problems varies locally in both 

space and time. The downstream boundary conditions can 
also create oscillations wl.en the mesh does not resolve 
the boundary layers. Correct treatment of the downstream 
boundaries eliminates this difficulty [95]. 

This provides an additional example of the ill effects 
of improper boundary treatment. As before the effects 
propagate into the interior and cannot always be easily 
traced back to its proper source. 
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3. HIGHER ORDER METHODS — FORMULATION 

As computer hardware improved, the need for more 
sophisticated numerical algorithms became obvious. In the 
early stages of software development a standard technique 
was a first order method while a high order method signi- 
fied a second order method. Computers were not sufficiently 
fast to consider two dimensional problems with fine meshes 
and so only low accuracy results were obtained. With the 
advent of faster and more structured computers it is now 
reasonable to achieve one percent accuracy for many two 
dimensional time dependent problems. Three dimensional 
problems are being solved with coarse meshes. With this 
situation it is necesse.y to analyze higher order 
methods. When one wishes accuracies of the order of one to 
five percent the higher order methods allow for a coarser 
grid than first or second order methods with no loss in 
accuracy. This coarser grid means that both computer 
storage and computer running time can be decreased without 
any deterioration in the solution. When even more accurate 
solutions are needed the advantages of the higher order 
methods are more pronounced. The results of [188] show that 
even with low accuracy requirements the fourth order method 
was more efficient. 

Higher order methods may not always be advantageous 
or feasible. Higher order methods usually require more 
computer time per time step than lower order methods. Hence, 
efficiency is increased only if a coarser mesh can be used. 
There are various circumstances where the mesh is 


12 


constrained by considerations other than accuracy and hence 
little is achieved by higher order methods. One case is 
when the geometry of the problem demands a large number of 
points. For example, if one wishes to describe the many 
perturbations on a real wing then one needs many more 
points than are needed for reasonable accuracy with a 
second order method. Another example occurs in meteorolog- 
ical flows over the globe. The accuracy of any algorithm 
is limited by uncertainties in the physics of the model and 
in observational data. However, one cannot choose too 
coarse a grid or the topography of the earth is distorted. 
Similar situations occur in other fields where the basic 
equations being integrated have only limited validity. 
However, for the majority of cases where the mesh is 
constructed mainly on accuracy considerations the use of 
higher order methods can lead to large savings in time and 
storage. Furthermore, the implementation of these methods 
frequently does not require large modifications to existing 
codes . 

The construction of higher order methods has proceeded 
along one of three lines. Either extensions of existing 
finite difference methods, or finite element methods or 
spectral methods. Multidimensional finite element methods 
have not proven very successful for hyperbolic problems. 
There are problems with the inversion of a large 
unstructured but sparse matrix and have thus far not 
competed successfully with standard alternating direction 




techniques or explicit methods. Furthermore, finite 
element methods do not automatically yield stable boundary 
treatment (see e.g. [72] and [80]) and so much of their 
justification for elliptic problems does not generalize. 

In this section we shall discuss finite difference methods 
and spectral methods . 

High order finite difference methods have the 
advantage that they are similar to lower order methods. 

Hence, their implementation is easier and usually does not 
require major modifications to existing programs. Spectral 
methods are even higher order and frequently are "infinite" 
order methods. They have limited applicability to problems 
with complex boundaries (see however [151]) and their suita- 
bility for shock problems needs further investigation. It is 
possible to construct unconditionally stable spectral methods; 
also boundary conditions are more straightforward with spec- 
tral methods. We shall concentrate on the so called pseudo- 
spectral or collocation methods as they have wide range of 
applicability. Galerkin spectral methods are more costly 
because of the need to calculate convolution sums. Hence, 
they are limited to equations with a quadratic nonlinearity 
in which case fast methods are available to calculate the 
convolution sums [149]. Even in this case they are two to 
tnree times slower tl.-” a pseudospectral method. 

We first consider the one-dimensional equation 




(3.1) 
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and some extensions of standard second order methods. 
The leapfrog approximation to (3.1) is given by 


(3.2) 


u 


.n+1 


3 





) 


with X = At/Ax. A fourth order extension suggested by 
Kreiss and Oliger [111] is 


n+1 n-1 X 

(3.3) u. = u . - g 




(f*^ - f*^ ) 1 

^ j+2 j-2'^ 


At the boundaries (3.3) can be supplemented by 
(see [64] and [143] 




X 

3 


(-llfQ+18fJ-9f2+2f3) 


(3.4) and 



n+1 n-i X ( -,n ^^n 

' “l '3 ^^1* ”2' ^3) 


n-1 


^ Xp r n+1 » n^ n-1 

2 ^'"l ■ "^1 > 


with p _> spectral radius of A = Sf/9u. 

These equations can be trivially solved for and 

Similar expressions can be derived for and 

To these conditions must be appended the appropriate 

g f 

boundary conditions. This scheme is stable if X 1 ‘72 
Another standard scheme is the Lax-Wendroff method. 

A two step version proposed by MacCormaclc [124] is 


(1) n , / £n ^n. 

“3 = "3 ' “'j+i ■ 

i '“3 " 




(3.5) 
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Another variant uses backward differences for the 
predictor and forward differences for the corrector. A 
fourth order extension is given by Gottlieb and Turkel 


[ 69l« [188]. 


u : 
3 


( 1 ) _ 


(3.6) 


u 


n+1 


u” . i (7f5 - 


T S ^ f 


The boundary treatment for the predictor is 


(3.7a) 


u 


(1) 

N-1 


- 1 '«[! 

- Ci- 


(1) 

N 






n 


4f" ) 

N-3 


and for the corrector 


u 


n+1 


2 6 


(3.7b) 


n+1 


u. 


(4f^^^-17f2^^+28f|^^-15fjj^^ ) ] 


1 . n. (1) X . .(1) . .,{!). ,(1) ..(Dx, 


The boundary conditions (3.7) are identical to using cubic 
extrapolation to calculate the fluxes exterior to the 
domain and then using (3.6). There exists another variant 
with backward differences in the predictor and forward 
differences for the corrector. For fourth order accuracy, 
in space, one must alternate the two variants. This scheme 
is stable if A < 2/3. 

d U — 


The schemes that have been presented are second order 
in time. It is possible to develop four step me the a s that 
are fourth order in both space and time. One such method 
is 


(1) 1 / n . n. 



n 

j+1 



(3.8) 


n , n ^ n 

u^2) = _2+l 2 2ll - A (fd) - f(l) ) 

3 8 2 3-r 


n 


n 


n, 


D + i 


-^V2-^Vi)-^^(yi^^? 

8 


t A [f 12) . £'2),l(.fn^2 4 3f - 3f £j.^) 1 


n+1 
u . 


= [f<^l-f(3Uf<2) -f(2) 

3 6 3+4 3-4 3+1 3-1 


+ |(-f + 7f ^^^1 f^^2^ 




1 +, 


^-2 




+ - + lOf" T - lOf" 1 + , 

16 j+2 ]+l ]-l j-2 


] . 


8 f 

This scheme is stable if A-;^ < 1 . 

9u — 

Other variants are given in [ 1 ] and extensions to 
multidimensions are described in [186]. Third order methods 
in time are given by Burstein and Mirin [30] and Rusanov [164], 
[165] . Steppeler developed third order methods in space 
and time based on an explicit evaluation of the Taylor 
Series combined with a third order finite element approxi- 
mation [175], [176]. 
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An alternative approach is to use Richardson extrapola- 
tion at each time step to increase the local accuracy in 
time. For the leapfrog method, which is symmetric about 
t +At/2 , this is straightforward while for the multistep 
methods it is more complicated. This approach increases 
the domain of dependence and also decreases the time step 
allowed by stability. It has the advantage of simplicity 
and of generalizing to muitidimensions ([187]). 

The schemes considered until now have been explicit 
methods that are fourth order in space. We now present 
higher order implicit methods. Comparisons between explicit 
and implicit schemes are presented in Section 5. We 
consider the compact implicit scheme [149], [202]. 

(3.9) (^)Au^_j^ + aAu^ + 

= - 1 - f^l)l 

with Au^ = - u^. 

: 3 3 

If a = 1 we have the standard implicit schemes. These 
are second order in time when C = 1/2 and unconditionally 
stable for C ^ 1/2. If a = 2/3 the schemes are fourth 
order in space. To solve (3.8) requires an iteration pro- 
cedure. An alternative is to expand Let A = , 

then a substitute for (3.9) is 


(3.10) (1^ -XCa")Au^_^ + aAuj + (^2. +ACA")AUj^j^ 


= - - (f - f ) 
2 '^j+1 j-l' 
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This scheme has the same formal order of accuracy in space 
and time as (3.9). Experience has shown that if one is 
approaching a steady state or else if the time evolution 
is slow then (3.10) is accurate even with large time steps 
(see [14]). However, for some problems the iteration proce- 
dure may be necessary. For sufficiently large time steps 
approaching a steady state and ^ = 1 it can be shown that 
(3.10) is approximately Newton's method for solving nonlinear 
equations [105]. 

Using a trapezoidal rule with end corrections it is 
possible to construct a fourth order in space and time 
formula that requires only three mesh points in the x 
direction. Let A be the Jacobian of f with respect 
to u. Then (3.9) can be replaced by 


(3.11) 


|(-4au5_^t Au^-4Au"^^) 


_A r^n+1 ,n+l 
4 ^^j + 1 “ ^j-1 


_ L rf“' - f 


A^r^n+l,,n+l r-n+1, .n+l,,n+l ^n+l. 


+ A.,i(f.,,-f.)-A. ,(f.-f. ,)] 
3+h 3+1 3 3-i 3 3-1 


If one wishes to linearize this scheme with fourth order 
accuracy in time, one must first calculate a second 
order accurate predictor. Hence, the generalization of 
(3.10) requires a predictor and corrector step. An 
additional difficulty with (3.11) is that the matrix 
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A ^ 

inversion is stable only when A~ ^ 1 . Hence# even though 
the scheme is uncon<?itionally stable one must restrict the 
time step in order to get a well conditioned matrix problem. 
Furthermore, it is difficult to generalize (3.11) to multi- 
dimensional problems since an alternating direction method 
will reduce the time accuracy of the algorithm. This scheme 
as well as some noncompact higher order implicit methods 
are described in [92]. Fourth order, in space, methods 
for equations with both hyperbolic and parabolic terms are 
considered in [41] and [161]. Some comparisons for the 
boundary layer equations are presented in [205]. It is 
also possible to construct higher order methods for these 
problems by combining high order Dufort-Frankel and leapfrog 
methods [68], [111]. 

Since these schemes are compact, special schemes are 
only required at the boundary. The author has found that 
the box scheme is especially accurate even though it 
reduces the order of accuracy. Hence, a boundary treatment 
for (3.10) would be 

(3.12) ( 1 - - ^ - ^ - *^ ..)Aug + (i±-|^)AuJ . -X(fJ-fJ) . 

Further discussion of implicit methods is presented in 
section 5. 

All the above methods are extensions of standard finite 
difference formulas. Hence, they have the advantage that 
it would not involve much programming effort to change an 
existing code. However, although they are of higher order 
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they are still only fourth order accurate in space. An 
alternative is to use a pseudo-spectral scheme. For suf- 
ficiently smooth solutions these methods converge faster 
than any power of the mesh size (see (70]). 

For periodic boundary conditions the most appropriate 
expansion series is the Fourier series. To solve (3.1) 
with a leapfrog-lilce method in time one algorithm is 

(3.13) u . = u . - 2 (Atf ) . 

J j * J 


where 

(3.14) 




I i£^G(kstl 
k=-N " 


•J 

where G (1^ At) = Ic At + 0((k At)"^). Given f this requires 
two fast Fourier transforms. A good choice for G is 

(3.15) Giklt) - (8 sin Oca At) - sin ( 2)co At ) /6o 

3 f 

with c 1.4 - (spectral radius of ■^)- This scheme is 
unconditionally stable. 

For problems requiring some dissipation one can 
replace the leapfrog method by a Runge-Kutta method. If 
one wishes to use a splitting method for multi-dimensional 
cases then in practice one is limited to second order in 
time methods. Unfortunately, the second order Runge-Kutta 
method is unconditionally unstable for the Fourier method. 
Hence, we consider a modified method similar to ( 3.13 - 3 . 14 ). 
The modified Euler scheme is given by 
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u • u” - 

( 3 . 16 ) 
with 

(Atf )!? » ? f,^G, 

^ J k— N ^ 

( 3 . 17 ) 

(Atf ) . - f LG2(kAt)e''^^’'/^ 

^ ^ k— N ^ ^ 

with Gj^(kAt) - kAt, G 2 (kat) = kAt for small kAt. 
Choosing Gj^(kAt) » GjtkAt) * kAt the stability condition 
is AN^At £ 8 with A » 3f/9u [70]. Since this is very 
restrictive we substitute this with 

Gj^(kAt) = (-e^^+ 8e^-7)/6o 

(3.18) 

G2(kAt) = (7-8e“^+e"^^)/67 

where z = ikoAt and r > 1.4c'(A). As before the scheme 
is unconditionally stable with these parameters. 

If we consider the simplified equation 

(3.19) Uj. + a (x) Ujj = 0 , 


the standard stability proof for the Fourier method is 
valid only if a(x) doesn't change sign (see [57]). The 
correct analysis of (3.19) when a(X(j) » 0 depends on 
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distinguishing between mesh^stability and time-stability. 
The Fourier method for (5.19) is stable for 0 < t < T 
as the number of modes increases. Hence, by the Lax- 
Equivalence theorem the approximation converges to the 
analytic solution. However, for a fixed number of modes 
the error may increase exponentially as time increases 
(when sCXq) “ 0). Hence, the Fourier method in its 
original form may not be useful when a(x) changes sign 
(74]. To overcome this difficulty we replace Gj^(kAt) by 
Gj^(kAt)p(k) and similarl- for G2 in (3.21). We choose 
the cutoff p(k) so that it is one for small k and goes 
to zero for tho highest modes. Majda et al. [131] and 
Kreiss and Oliger (113] have shown that this modification 
to the Fourier method is stable. In fact, using these 
cutoffs with a C''ebysheY method (to be described) one can 
solve the Riemann problem of fluid dynamics and resolve 
the shock and contact discontinuity within one mesh width 
(D. Gotf lieb, private conimunication) . 

For bounded domains, ?n expansion in Chebyshev poly- 
nomials is more appropriate. Since the eigenvalues of the 
operator have a negative veal part, we consider a two step 
scheme rather than n :.-apfrog method. Hence, we have 


= u" - A (Itf )" 
X 3 


n+1 n 
u = u . - 




(3.20) 
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Let Xj ■ cos iTj/N, one can then expand 

N 

(3.21) “ I j-0,...,N 

J k«0 


where 170] 


^ ^ 'j ^ • =0 ■ =N ’ 2 

Cj^ « 1 t * 1, . . .N 


Equation (3.22) is solved using the FFT. We then have 


N 

(3.23) ). - I b.T. (X.) 

* ^ k»0 ^ 


with 


(3.24) \ ‘ I apG(£At) 

^ £=k-H 

£+k odd 

and G(£At) « £At + 0((At)^). In this case a reasonable 

■" Ol 2, A t 

choice is (see [73l)» z * e , and 


(3.25) 


G(£ At ) 


(.*. - 18z + 


9z^ - 2z^)/6a 


with a approximately equal to the spectral radius of 
The extension to multidimensional problems in 
rectangular domains is straightforward- The leapfrog 
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methods generalize in the obvious manner. For the Hunge- 
Kutta methods splitting techniques are avaiable for multi- 
dimensions while for the implicit methods alternating 
direction techniques are the most appropriate. 

We have stressed methods that are second order in time. 
There are several reasons why higher order methods in time 
are less appropriate. First, the extension to several 
dimensions is complicated as splitting and alternating 
direction methods are second order in time unless complicated 
algorithms are constructed. The necessity for a full multi- 
dimensional scheme usually results in a method which requires 
many operations per time step. For problems v‘’.ere the time 
dependent equations are being used as a relaxation scheme 
to a steady state there is no obvious advantage to a higher 
order method in time. Even for time dependent problems the 
important physical phenomena frequently move at speeds 
considerably slower than the fastest signal speed allowed by 
the equations. Hence, these motions are usually easy to resolve 
and only the variation in space is rapid requiring higher order 
techniques. The other major argument for second order schemes 
in time is that one can always increase the time accuracy by 
taking smaller time steps. The computer work involved varies 
linearly with the time step and so this may be more efficient 
than a costly high order method. Furthermore, choosing 
smaller time steps does not increase th*- storage. ilowever. 


for explicit methods doubling the mesh points increases 
the work by 2^^^ in d dimensions and the storage 
requirements increase by 2*^. Hence, the use of high 
order methods in space to limit the number of mesh points 
can be very advantageous. Furthermore, one can increase 
the accuracy in time by locally using Richardson extra- 
polation il87J. In general it is not known when higher 
order methods in space or time are more efficient than 
lower order methods. One case where the author has found 
higher order methods, in time, useful is instability 
studies. Time inaccuracies can introduce errors which 
cause artificial instabilities. As more experience is 
gained better guidelines should be available. 
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4 . high order methods — RESULTS 


In [150] and [186] a range of problems are considered 
to demonstrate the advantage of high order methods . Here 
we consider two of these cases which illustrate many of 
the phenomena of more complex situations. 

The first problem is flow in a nozzle. The equations 
are one dimensional but many two dimensional effects are 
included by specifying A(x) the area of the nozzle at 
position X. The equations of motion are 


(Ap)^ + (pAu)^ = 0 

(4.1) (Apu)^ + [ A(pu^+ p)]^ = A^p 

(AE)^ + [Au(E + p)]^ = 0 

with p = (y-l) (E - j pu^) and y = 1.4. 


The solution to the steady state equations is known 
for both smooth and shocked profiles [43] . Hence, we march 
the equations toward a steady state and we can then compute 
the errors between the computational solution and the known 
analytical solution. 

We linearize (4.1) about a constant state (Pq,Uq,Eq) 
and drop lower order terms. One then finds that the 
characteristic values A and eigenfunctions v are given 
by 


= u 
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(4.2a) 
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(4.2b) ^2 - “o ^^2 “ 


(4.2c) - Uq - Cq - Uq 
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pu + E 


u + E 


with 


Cq = y(Y“1) 


^0 ,2 


At the inlet the flow is subsonic and we specify 
p and E equal to the known steady solution. This is 
supplemented by a difference equation for the characteristic 
variable which is coming into the boundary. Given p, 

E and v^ one can trivially solve for pu. With the 
explicit method it was not found very important to solve 
for v^ rather than pu. However, with the implicit method 
and large time steps the use of (4.2c) was crucial. As 
pointed out before the use of noncharacteristic boundary 
treatment frequently manifested itself as a negative pressure 
at the sonic point. The only way of identifying it as a 
boundary difficulty was to observe that the difficulty dis- 
appeared when one solved for v^ rather than pu by the 
boundary difference scheme. 

For smooth solutions no boundary conditions are speci- 
fied at the exit. Hence, the special boundary treatment 
is used for all three variables. For a shocked profile the 
pressure is specified at the exit. Again, with the explicit 
schemes it was not found necessary to use characteristic 
variables. Instead we calculated o and pu by finite 
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difference approximations to (4.1). E was then calculated 
knowing p, pu and p. With the implicit methods one must 
first linearize the given boundary condition so that itera- 
tion is not required. This is similar to the manner in 
which (3.10) was derived. Then p = Pexit replaced 
by 

u^ _ n 

(4.3) ^ A p - UgAm + AE = ; A o = p^'^^-p" 

By inspection one verifies that in the steady state 
p = (4.3) is supplemented by two equations based on 

the box scheme (3.13) for the characteristic variables v^^ 
and V 2 * For the higher order methods A' (x) must be also 
calculated to the same order. Hence, for all problems both 
A(x) and A' (x) were given analytically. For the examples 
considered A(x) was chosen as a hyperbolic cosine. 

For smooth profiles, it was found that the major 
difficulty occurred at the throat where the flow is sonic. 
When an artificial viscosity is not used many of the 
methods gave rise to an expansion shock. For most of the 
problems the solution is considered at steady state if p 
changed by less than £ = 10 ^ in one time step. For some 
of the runs with high resolution a smaller t was chosen. 

In Table 4.1 we present the results for a smooth 
profile. The explicit MacCormack methods did not require 
an artificial viscosity while the implicit method did. 

When needed the artificial viscosity is added explicitly 


at time t even for the implicit methods. We see from the 
table that the fourth order explicit method is more than 
twice as efficient as the corresponding second order method. 
It also requires less than half the storage to achieve 
three digit accuracy in the steady state. For four digit 
accuracy the efficiency factor increases to 16. 

For implicit methods the accuracy and stability depends 
very much on the boundary treatment. In [72] it is proved 
that one does not need to include the given boundary condi- 
tions in the implicit method. Instead# after each time step 
one can correct for the boundary conditions. For equations 
(4.1) it was found that this worked only for time steps 
less than three times the Courant limit, otherwise 
nonlinear instabilities arose. For larger time steps it 
was necessary to incorporate the given boundary conditions 
in the matrix to be inverted. With time steps of about 
15 times the Courant limit steady state was reached very 
rapidly. In this case the higher order method was only 
about thirty percent more efficient. We speculate that in 
going to a steady state the use of the lower order box 
scheme at the boundary deteriorates the accuracy. For a 
wave equation the fourth order implicit method was 3-4 
times more efficient than the second order implicit 
method (see also [63]). For both these cases there is no 
improvement in either accuracy or in the rate of convergence 
when a second order in time method is used rather than a 
first order in time method. 
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This problem, (4.1), was also solved with the Chebyshev 
collocation method (3.16) - (3.17) with Gj^(kAt) * G 2 (kAt) « 
kAt. This is stable with time steps about .1 times the 
Courant limit. The error at the steady state is concentrated 
at the sonic point. A simple automatic postprocessor at 
the conclusion of the code removes this error. The result 
is an error level 20 times smaller than with the fourth 
order method, as seen in table 4.1b. 

We next consider (4.1) with a shocked profile. The 

solution to most of the methods is not a monotone function 

of the mesh. Hence, instead of comparing schemes for a 

given error tolerance we find the asymptotic rate of 

accuracy. This is found by a least square estimate based 

on about 50 runs with different meshes. If one includes 

the shock area in the error analysis then all the methods 

1/2 

converge like (Ax) in the L2 norm and as Ax in the 
norm. When we exclude a fixed physical distance about the 
shock the scheme behaves statistically according to the 
formal accuracy of the method. With the implicit methods 
instabilities appeared when time steps larger than five 
times the Courant limit were used. 

As the second example we consider a two dimensional 
problem in dynamic acoustics. The equations are given by 
(7.1) where the radiation boundary treatment for this 
problem is discussed. The numerical algorithm used to 
solve this problem is the second and fourth order two 
step schemes (3.5) and (3.6) together with a splitting of 
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the two dimensional problem into a series of one dimension- 
al problems. This code was run on the CDC-STAR at NASA 
Langley. More details of the algorithm are found 
in [10]. 

Two cases are considered here. In the first case the 
mean flow is zero and the forcing function is a delta 
function in space and harmonic in time; for this case the 
analytic solution is known. The fourth order method is 
more accurate than a second order method with twice as many 
mesh points. The efficiency gain is about 400 percent at 
large error tolerances and increases at lower error toler- 
ances . 

In the second case we use a realistic mean flow modeled 
after experiments and a harmonic source located two jet 
diameters downstream of the jet exit. Varying the mesh 
and boundaries demonstrates that the fourth order method 
with 12000 points yields essentially the analytic result. 

We then measure the peak acoustic pressure as a function of 
the angle. The fourth order method with 8800 points gives 
much better accuracy than the second order method with 
16000 mesh points. The fourth order method is now being 
routinely used to solve problems with a variety of sources 
including pulses, convecting sources and quadrapoles. Second 
order methods would not be feasible for these problems because 
of storage and cost limitations. 

In general we have found the fourth order methods to 
be three to five times faster while still giving accuracy 
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conparable with the second order method. The savings in 
storage is a factor of two per space dimensions. For 
error levels less than one percent the efficiency of the 
higher order methods increase. Similar conclusions were 
reached in [150] . Coding changes for the finite difference 
methods are minimal given a second order program with 
either an explicit or implicit scheme. 

On some simple problems spectral methods behave even 
more efficiently than the fourth order schemes. However, 
they require coding a new program with careful attention 
being paid to the implementation. At present there has 
been little experience with these methods for large scale 
hyperbolic problems especially in complex geometries or 
with nonsmooth profiles. 
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5. IMPLICIT METHODS 

The use of implicit methods to solve hyperbolic 
equations has been increasing in recent years (see e.g. 

[20]f [87])« The rationale behind this is that implicit 
methods are frequently unconditionally stable. Choosing 
large time steps can then more than compensate for the 
extra work per time step. 

Even though there are no stability restrictions on 
the time step nevertheless the time step is still restricted 
by accuracy requirements. We consider the simplest equation 


(5.1) “t " ° 0 < X £ 2tt 


which we approximate by the Crank-Nicolson scheme 


(5.2) 


n+1 


- V . = 


At 

4Ax 




Assume Vj = sin(jkAx). Since the Crank-Nicolson formula 
(5.2) is nondissipative the only errors are phase errors, 
so that = sin(jkAx+a) . For the differential equations 

a. = kAt. The numerical phase for (5.2) is given by 


v5.3) 


‘N 


= sin 


-1 


sin(kAx) 


1 + (At) sin (kAx)/4(Ax)^ 


In Table 5.1 we present the phase errors 

kAx = Yo> f* phase errors increase 

dramatically as we choose At/Ax much larger than one. 


35 


This is because the Courant condition dt/Ax £ 1 has 
physical significance beyond a fomal stability condition. 
With larger time steps, we are not following the wave 
correctly. 

The justification for implicit methods arises only 
by considering equations more complicated than a wave 
equation. A simple model is provided by 

(5.4) *^t '*x * A(l-p)ccs (x-pt) 

A solution to (5.4) is 


(5.5) u(x,t) ■ A cos(x-pt) + B cos(x-t) 

If p << 1 and B << 1 then accuracy requirements only 
demand 


(5.6) 


Ax 


1 


which is much weaker than the stability criterion for 
explicit schemes At/ Ax 1. 

Another example is provided by systems of equations 
with widely separated eigenvalues. For simplicity we con- 
sider the uncoupled system 


u 

V 


t 

t 


+ au. 


+ bv. 


0 


(5.7) 


0 
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with u • A sin(x-at), v ■> B sin(x-bt) and A << B, 
a >> b. As before explicit achemea require aAt/Ax ^ 1. 
However, aince A << B accurate aolutiona only require 
bAt/Ax <, 1. 

Both of theae exa»:tples demonatrate the phenomena of 
different time scales. For these problems the tin^ step 
of an explicit method would be limited by the speed of 
the fastest possible mode. For implicit methods the time 
step is chosen to resolve the slower modes which carry 
most of the energy. Both these simple illustrations have 
many practical applications. For example a slowly oscil- 
lating wing or rotor will induce wave motion with much 
slower speeds than that of free motion. In meteorology 
or plasma physics the usual speeds of propagation are 
much smaller than the fastest signal speed. These will be 
discussed further in chapter 8. The possible extensions 
of explicit methods for these problems will also be dis- 
cussed. 

A similar justification for implicit method occurs 
when the time step for an explicit method would vary 
dramatically between different regions. An extreme example 
occurs in laser fusion where the diffusitivity can change 
by many orders of magnitude across the domain. In magnetic 
fusion the existence of near vacuum regions create areas of 
very high speeds compared with the center of the 
plasma. In other problems similar effects occur due to 
different mesh sizes in different regions. For problems 
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in Mtarology the time step for an explicit echene is 
determined by the email mesh voIukms near the pole due 
to converging latitude lines. In Navler-Stokes problems 
with boundary layers extresMly small iMsh sizes occur 
nMr the body In order to resolve the boundary layer. 

In all these problffias the time step of an explicit scheme 
is governed by a small region which may not be the area 
requiring the greatest resolution in time. 

A second type of problem suitable for Implicit 
method are those cases for whicn only the steady-state 
solution is desired. The time dependent equations are 
used nwrely as a device for obtaining iterative solutions 
to the steady state equations. In this case there is no 
need for the numerical method to accurately follow the 
transient. Indeed, one way to accelerate the iteration 
process is to make the scheme inconsistent with the transi- 
ent solution. We only need guarantee that the numerical 
steady state achieved is independent of the iteration 
procedure. Consider the general equation 

(5.8) u^ ■ Lu . 

A simple way to ensure the correct steady state is to 
solve for u” ■ u”^^-u” at each time step. The algorithm 
then has the form 


q”au” 


Lu‘ 


(5.9) 
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The steady state is defined by - 0 and so we have 

Lu" *0. Since the time accuracy is not important we 
do not wish to iterate at each time step. In the interest 
of efficiency q" should not depend on the unknown variable 
In many cases this can be achieved without loss of 
accuracy by a linearization procedure such as introduced in 
[14], [23], and [122]. Higher order (in space) of these 
schemes were defined by (3.10) and (3.11). 

If we use the backward Euler formula for (5.8) we have 

(5.10) • AtLu”'*'^ . 

Linearizing f this equation we have that 

(5.11) ^d-Atj")Au" - Lu" , 

J is the Jacobian of L with respect to u. As At 
increases this scheme approaches the Newton-Raphson iteration 

(5.12) -j”Au" - Lu" 

This is a property only of the backward Euler formula and 
is not true for other selections of parameters in (3.9). 

In general it is not necessary that o” be a function 
of the exact Jacobian of L since we are only interested in 
the steady state. This freedom can be utilized in several 
different ways. One possibility is to change in such a 
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way as to simplify the solution of the linear set of 
equations (5.11). In most cases involves block 

tridiagonal matrices. If these raatricec can be inverted 
by Gaussian eliminations without pivoting then the in- 
versions can be accomplished by the Thomas algorithm in 
O(m^N) operations where m is the block size and N 
is the number of mesh points (see 11031). Hence, for 
general r. ’ jck tridiagonal matrices and m larger than 
3 or 4, most of the work is in inverting the full blocks. 
When the fluid dynamic equations are written in velocity 
form these blocks can be decomposed as a direct sum of 
smaller blocks and so the process can be speeded up [23] . 
When the momentum form of the fluid equations are used 
full blocks occur. However, by using the freedom in q” 
one can simplify these blocks at the expense of making the 
numerical scheme inconsistent with the time dependent 
equations [174]. In fact even when solving the steady 
equations directly using Newton's method [17] J need 
not be the Jacobian of L. This generalization leads to 
the use of quasi-Newton methods (for a survey see [46]). 

Alternatively one can choose in such a manner as 

to speed up the convergence. McDonald and Briley [128] 
consider methods with different It at different mesh 
points. This can be viewed as a matrix conditioning of 
the linear equations. This is especially promising for 
parabolic equations in two space dimensions where the theory 
of parameter selection for A.D.I. methods is well developed. 
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Implicit methods have the disadvantage that they require 
the solution of a large number of coupled equations at each 
time step. Therefore, the reduction in the number of time 
steps compared with explicit methods must be weighed by 
the increase in the number of arithmetic operations 
required for each step. To simplify the inversion process 
in several dimensions alternating direction methods are 
generally used. In this case one needs to invert a block 
tridiagonal matrix for each direction. As mentioned pre- 
viously this is very expensive when the block sizes get 
large. For example, in the magnetohydrodynamic case the 
blocks are 8x 8 matrixes. In combustion problems there 
is at least one partial differential equation per species. 
Hence, for complicated chemical processes very large blocks 
can be generated. These situations render standard implicit 
methods impractical as the work increases with the cube of 
the block size. In some cases one can use knowledge of the 
block structure to reduce this work [23]. However, when 
shocks appear and the conservative forms of the equations 
must be used, full matrices are unavoidable. 

Furthermore, the Thomas algorithm is an inherently 
serial algorithm and so inefficient for many vector pro- 
cessors. In three space dimensions one can perform many 
tridiagonal inversion simultaneously to partially vectorize 
the procedure. However, this demands large storage require- 
ments. The use of a cyclic reduction method is more efficient 
for a vector machine but is still very far from optimal for 
these pipeline machines. Hence, many of the advantages of 
A.D.I. are negated on machines as the STAR-100. 
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An additional difficulty with alternating direction 

2 

methods is that they introduce 0((At) ) perturbations 
into the matrix q" of (5.10). In this case the scheme 
no longer approximates the Newton method for large At 
and consequently there is a reduction in the rate of con- 
vergence to a steady state. Marching to a steady state 
using large time steps one wants to use the delta form 
(5.10) so as to ensure that the steady state is independent 
of At. In two space dimensions the alternating direction 
methods which solve for u^^^ and Au*^ are equivalent, 
but in three dimensions they are not. The three dimensional 
algorithm is unconditionally stable in the linear case if 
one solves for u but the steady state depends on At. 

On the other hand if one solves for to produce a 

steady solution independent of At, the algorithm is uncon- 
ditionally unstable for scalar problems. As the entropy 
equation is essentially a scalar equation this method has 
difficulties for many systems. Only the addition of 
viscous terms can stabilize the procedure. 

* 

Several alternatives have recently been advanced as 
substitutes to alternating direction methods. Steger and 
Warming [176] have suggested splitting the flux vector into 
two parts corresponding to the positive and negative eigen- 
values. Each part is then solved using the one sided 
differences appropriate for the corresponding eigenvalue. 
Jameson and Turkel [105] have proposed a method based on a 
LU decomposition. In this method the lower and upper factors 
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are chosen for computational ease. The resulting scheme 
can be chosen to be second order accurate in space. It is 
demonstrated that the crucial ingredient is that each of 
the L and U factors be diagonally dominant. The dia- 
gonal dominance of the final scheme is irrelevant. This 
scheme is stable for one, two and three dimensions. A 
common feature of both these schemes is that only two factors 
appear independent of the number of space dimensions. For 

one dimension this is a disadvantage since it introduces 

2 

perturbations of order (At) and so '.l^ws dov^n the con- 
vergence rate. However, for three space dimensions the 
A.D.I. schemes have changes of order (At) ^ from the baclc- 
ward Euler method and so converge to a steady state even 
slower than these methods for large At. The requirement 
of three sweeps through the mesh for a three dimensional 
A.D.I. method is also a disadvantage when not all the infor- 
mation can be stored in core. 

Since the bacl^ward Euler method is a good approximation 
to Newton's method it may be advantageous to use this method 
even for multidimensional problems. The resulting matrix 
is no longer tridiagonal and hence it is necessary to find 
some efficient method to invert the matrix. Band Gaussian 
elimination solvers require excessive core especially since 
pivoting must be used. For parabolic problems that arise 
in laser fusion Kershaw [108] has used a conjugate gradient 
method for inverting the matrix. Similarly Orszag [151] has 
advocated the use of conjugate gradient to invert the full 
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matrices that arise from spectral methods. Brackbill [20] 
has used a SOR method in conjunction with his implicit 
method for the nonlinear MHO equations. 

Until now we have concentrated on the rate of con- 
vergence to a steady state. An equally important topic 
is the accuracy in the steady state. We wish to stop the 
iteration process in such a manner so that the error in 
the iteration process is below the truncation error of the 
steady state scheme. In using (5.10) we have assumed that 
in the steady state Au’’ =0. In many codes, one iterates 
until Au^ is below some given error tolerance, e. We 
then have that Lu*^ = Q^e so that the error in Lu is 
effected by the operator, Q. Therefore, it is better to 
use the residual Lu*^ as a measure of the steady state 
rather than Au^. Since, L is in general not an elliptic 
operator it is difficult to measure the deviation of u*^ 
from the steady state, even given that Lu” = e. It is 
also important to choose initial conditions that are a 
reasonable guess to the steady state solution. 

In summary, implicit methods have been successful when 
one is careful to match the physics with the method. These 
methods are less appropriate for wave-like equations where 
one wishes to follow all the possible modes of propagation. 

Attention to boundary treatment is even more important 
for implicit methods than for explicit methods. This is 
mainly because one wi.shes to use the implicit methods with 


large time steps. In many calculations one could use any 
one sided difference to some equation to supplement the given 
boundary conditions as long as small time steps were used. 
However, for large time steps it was essential, for 
nonlinear stability, that the characteristics variables, 
e.g. (4.2), be used. Skollermo [172] gives examples where 
the use of an explicit boundary method can force a stabi- 
lity condition for the entire method. 

For complicated flows the choice of a poor boundary 
treatment may be difficult to judge. For example, in- 
correct treatment of outflow boundaries will severly slow 
the rate of convergence to a steady state. Rudy and 
Strikwerda [162], [163] and Thomas [182] demonstrate 
that overspecification can be particularly inefficient. 
However, if one judges the results by comparison with 
experiment one would never sense the incorrect boundary 
treatment. Gustaffson and Kreiss [85] show that for this 
case the steady state may depend on the initial conditions. 

In any case it would not be obvious that the slow rate of 
convergence is due to the boundary treatment. Thomas [182] 
describes boundary treatments for other types of boundaries 
that occur in Navier-Stokes flow (see also [4]). 

Due to the unconditional stability of implicit 
schemes it is not clear how to choose the time step. One 
simple procedure is 


(At) 


n+1 


(At) 


l|Au» II 

l|u^ II 
l|Au"-^H 
l|u"-^|| 


4S 


Au 


n 


u" - u”-l 


for some choice of norm. This has been used successfully 
in plasma diffusion problems (D. Nelson, private communica- 
tion) . A more sophisticated choice is to compare t%iro 
iterates for time t as is done in o.d.e. solvers. When 
going to a steady state At can be viewed as an interation 
parameter, as discussed previously. 


TABLE 5.1 

Phase errors for the Crank-Nicolson method 
as a function of the time step. 
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6. MOVING BOUNDARIES AND ADAPTIVE GRIDS 

In the usual fluid dynamical situation boundaries are 
considered as fixed in time. In many situations, especi- 
ally for analytical results, the fluid is considered as 
being confined to a rectangular or circular region. However, 
in many circumstances one must include the motion of the 
boundary as an important element of the problem. The 
movement of the boundary arises from many different 
factors which require different methods. 

The simplest situation arises when the boundary moves in 
response to an external force. This may represent a moving 
piston, a diaphragm or similar devices. The next 
situation occurs when the boundary represents a free 
surface. In this case the boundary represents the separa- 
tion between the domain of interest and some other region, 
for instance a vacuum or the general atmosphere. The outer 
region presents no resistance and the boundary moves in 
accordance with the forces exerted on it by the interior 
material. This occurs when metals are subject to a high 
temperature or pressure and begin to flow. Other examples 
occur when liquids are not in a container, as in 
water over a dam or water waves or raindrops. The most 
difficult problem arises when the moving boundary represents 
an interface between different materials. In many cases 
these are materials subject to the same differential equa- 
tions. The two materials differ only in their density or 
other material properties. The simplest such case is a 


contact discontinuity in fluid mechanics where both sides 
are the same material but jumps occur across the moving 
surface. In other situations the media on the two sides of 
the moving boundary are represented by different 
differential equations. Such conditions aris-; when an 
explosive gas impinges on a solid material. Another 
example is the interface between a plasma and a vacuum. 

In the latter case hyperbolic time dependent equations 
describe the motion of the plasma while the magnetic field 
in the vacuum is given by a time independent Poisson equa- 
tion. In these cases the boundary moves as a result of 
imbalance of forces from the two sides. 

The standard techniques to solve such problems are 
Lagrangian methods (see e.g. [19], [20], [96], [197]). 

With such schemes the boundaries are coordinate lines; 
this simplifies the algorithm. However, Lagrangian methods 
have several drawbacks. They are usually low order methods 
especially in regions where the mesh is nonuniform. When the 
motions are large the grid undergoes severe distortions which 
require rezoning the mesh. This rezoning is quite difficult 
in three dimensions. The rezoning usually results in a loss 
of mass and so should not be done too frequently. In 
addition, the rezoning formulas are usually not automated 
and require intervention by the user. The main drawback of 
Lagrangian methods is that they are very complicated and 
not user oriented. One method of simplifying the rczoning 
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difficulty is to use triangular meshes. At each time 
step reconnections are made when appropriate. This 
method developed by Fritts and Boris [60] has been 
applied to problems involving high shear. As a dual 
to this method is a scheme devised by Peskin [153] which 
is a grid free Lagrangian method. This latter scheme 
has been applied to the incompressible Navier-Stokes 
equations. A difficulty with both these methods are 
that they are difficult to couple with implicit time 
algo;..ithms and also the extension to three dimensions 
is computationally complicated due to the many possible 
configurations . 

At the opposite end one can use a strictly Eulerian 
approach and integrate across the boundary. To prevent 
smearing of the interface one adds some artificial compres- 
sion after each time step [91] . The location of each 
material is identified by a color function which is 1 in 
one region and 0 in the other region. This color function 
satisfies a convection equation which itself must be 
solved numerically without smearing. This method has not 
been used extensively on large scale problems and its 
applicability for multidimensional problems is question- 
able. It would be difficult to implement at interfaces 
where jump conditions need to be satisfied. These jump 
conditions depend on the physics of the situation and 
cannot be derived just from properties of the differential 
equation . 
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One alternative is to map each region onto a rectangle. 
If this mapping is done at each time step then one must 
interpolate between the grids at successive time steps. 
Instead, one can initially find a grid by any package, 
e.g. 1184]. This gives new coordinates C “ C(x,y), 
n “ n(x,y). We now allow the new coordinates to vary 
with time, so that C “ C(x,y,t) and n n(x,y,t). Given 
the differential equation 

(6.1) ''t ^x ° • 

This gets transformed into (see [154], [193]) 

( 6 . 2 ) = 0 

with 


(6.3a) 

and 

(6.3b) 


q » I W 

P = I(C^.w+ Cjjf + t^g) 
G = I (n^w+ Hyf + n^g) 

^t ~ ~ ^t^y 


I - x,y„ - 


CyHx 


(6.3c) 


j = 


l/l 
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In order to determine all these quantities we require 
the grid point speeds . One method is to construct 

a new grid at the advanced time. This can be done either as 
a function of gradients in the problem [50] or by solving an 
elliptic equation. Given (x,y) at t and t + At (x^»y^) 
can be calculated. An iterative procedure would be more 
costly but also more stable. Hindeman et al. [98] prefer 
differentiating the elliptic equation for (x,y) with 
respect to t. This yields a linear elliptic equation for 
(x^,y^) . Given (x^,y^) on the boundary this equation is 
solved at each time level. It is also important to solve 
for the gird in a manner which is consistent with the 
numerical solution of the differential equations [183] . 

These procedures require the solution of an elliptic 
equation at each time step. It is not clear that the over- 
head required by the mapping justifies its use. In addition, 
we would like to use information about the gradients of 
the solution to construct the grid at the new time level. 

At this time it is not known what are reasonable ways of 
accomplishing this especially for multidimensions. For 
example, it is well known that we can not allow the image 
of a square in (C,n) to become too distorted or else in- 
accuracies and instabilities may occur. Also, if the grid 
changes too rapidly in time one would expect difficulties. 

An experimental program for a parabolic problem is presented 
in 150] while one for hyperbolic systems is presented in 198). 
Oliger (146) and Yanenko (207) also investigate adaptive 
grids from a more theoretical viewpoint. 
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The differential equations are solved on the fixed 
Euler ian TOSh. At mesh points far from the boundary the 
standard Eulerian schenw is used. At grid points near the 
boundary one sided differences are used. Alternatively the 
fluxes can be extrapolated to the outside of each domain 
and then the standard scheme is used. The boundary itself 
is identified by a series of Lagrangian tracer particles 
which are allowed to move through the fixed Eulerian mesh. 
This boundary is used only to keep the regions separate and 
prevent diffusion of one material into another regioii. The 
only communication between different regions is via jump 
conditions across an interface. No differential equations 
are integrated across a boundary. 

As an example of the difficulties encountered with 
moving surfaces we discuss the impact of materials at high 
speeds. Metals impacted by gases or by other metals 
are subjected to high pressures which will deform the 
metal. The metals display elastic-plastic deformations. 

The differential system for these situations is given by 
the Prandtl-Reuss equations [97]. As usual, p, u, v, e 
represent the density, velocities and internal energy 
respectively. The total stress are given by 

are stress deviatories and p is the thermodynamic 
pressure. When the deviatories are less than the stress 
limit, i.e. ^ tc the flow is elastic. For 

A J 

2 2 2 
[st j = < the flow is plastic. < is constant in the 

simplest of models but is a function of various dependent 
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variables when work hardening is included. We shall con~ 

aider these equations wit:', cylindrical symmetry. 

|1 ^ 

Let - + u + V -jp . Then the equations in 

the generalized elastic regime are given by 




du 


^ du . ap 

Pat H 


at dr 


^ 1^1 
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with the strain rates 


ij 


given by 
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11 


if- ^ 3v vl 
" 3z " 3r " rj 


(6.5) 


•12 


1 9u . ^ 

2 dr 3z 


_ 1 f- ^ iH Y. 

^22 “ I 3r “ 3z " r 


To this we append an equation of state 


( 6 . 6 ) 


p = p(p,€ ) . 


Since these equations are hyperbolic we can con- 
struct characteristic equations for which differentia- 
tions in only two space-time directions occur. In con- 
trast to the fluid dynamic equation there are two speeds 
of propagation. The first is the shear or transverse 
speed given by 

2 

(6.7a) Cg = y/p . 

There is also the compressive or longitudinal speed 
(6.7b) c^ = c^ + 4 m/3p 


where 
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(6.7c) 

dp 

is the fluid dynamic sound speed. In addition to these 
two conoids there are three quantiti'^s that propagate 
along stream lines. Hence, in total there are seven 
sets of bicharacteristics. 

In the interior we solve (6.4) by any dissipative 
conservative scheme. We do not follow the shocks explicit- 
ly but instead capture them. The elimination of shock 
fitting eliminates many of the complexities of the problem 
and allows us to concentrate on fitting the interfaces. 

Since shocks are compressive it is possible to follow them 
without shock fitting. However, contact discontinuities 
will be unacceptably smeared unless some special procedure 
is followed. Hence, all free surfaces or interfaces are 
followed explicitly to prevent their diffusion. 

Due to the complications of the boundary we have 
chosen a simple scheme for the Eulerian mesh. All the 
results were obtained using a two step scheme devised by 
Burstein [29] This scheme uses data at a nine point 
rectangular lattice at the time level t to advance to t+At . 

At points far from the boundary the computation is straight- 
forward. At mesh points near the boundary the fluxes at arti- 
ficial points outside the domain are found by extrapolation. 
Quadratic extrapolation including the values at the 
boundary is used because of the complexities that arise 
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with an arbitrarily shaped boundary. After the flux 
extrapolation the nine point scheme can be used without 
any complications. Each side of an interface is updated 
independently. In no case are values on one side of the 
boundary used to update variables on the other side. In 
many cases different grids are used in different regions 
so that it would be difficult to integrate across an 
interface even for quantities that are continuous across 
the interface. 

All boundaries, which move through the fixed Eulerian 
mesh, are marked by Lagrangian tracer particles with 
position (z,r) and velocity (Ug.v^) . For interfaces between 
two regions there are marker particles on both sides of the 
moving boundary. The motion of these particles is governed 
by the equations 


( 6 . 8 ) 


dt ~ Ug(z,r,t) 
if = Vg(z,r,t) . 


Tie velocity (u„,v„) is the local fluid velocity and is 

DO 

found by extrapolation from the interior. Across an inter- 
face only the normal velocity is continuous. Hence, there 
will be a tangential slippage of the position of the tracer 
particles on one side of the interface with respect to the 
other side. The system of equations (6.8) are solved for 
the new positions z(t * i.t), r(t + .’.t) at each boundary 
point by a first order method. 
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Since each tracer particle is moved independently the 
distance between these particles will vary. In many cases 
the particles will bunch up in some regions with very 
few boundary particles in other sections. To prevent this 
the particles are rezoned if the distance between two 
particles differs by more than 25 percent from the average 
distance between particles. This rezoning consists of 
creating a new set of tracer particles which are equally 
spaced. Any information needed at the tracer particles 
are found by interpolation. All distances and interpola- 
tion formulas are calculated in terms of arc length along 
the boundary. This is a one dimensional operation and so 
much simpler than a full Lagrangian rezone of the entire 
two dimensional grid. Furthermore, since the interior 
values are unaffected by this transformation, the rezoning 
procedure cannot affect the conservation of mass, momentum 
or energy. 

In order to advance the solution we need the dependent 
variables at the boundary itself. It is here that the jump 
conditions affect the solution. We first describe the 
boundary conditions at an interface separating two contigu- 
ous elastic domains. The physical laws that apply at 
material interface boundaries are (1) continuity of the 
normal velocity, u^ , (2) continuity of the normal stress, 

, and (3, 4) the specification of the shear stress. 


^ng ' 


on each side of the interface as a given function of 
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the normal stress. Here we only consider the condition 

that T = 0 on each side of the interface, which is the 
ng 

free slip condition. One could also consider welded sur- 
faces where laws (3, 4) are replaced by the condition that 
the entire stress matrix be continuous across the inter- 
face . 

It may be verified that, for the elastic equations, 
there are four characteristic waves emanating from the mov- 
ing interface, i.e., from each side of the interface two 
waves propagate away from the boundary. Hence, at the 
interface boundary, a total of four conditions need to be 
specified, so that by satisfying the above physical laws at 
the interface boundary, the boundary motion can be determined. 

To implement the boundary conditions we consider a 
coordinate system with coordinates, n the normal and s 
the tangent to the boundary at each marker point along the 
boundary. Let S(z,r) be the deviatoric stress tensor as 
a function of (z,r). Let R be the rotation matrix 


(6.9a) 


R 


cos 0 
-sin 6 


sin 6 
cos e 


with 0 the angle measured clockwise from the Z axis 
to the normal. Then 


S(n,s) = R S (z,r)R^ 


(6.9b) 
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is the rotated deviatoric tensor. In this reference 
frame, the stress tensor T with components can 

be computed. The velocity vector is given by 


(6.9c) 


’ • 

u 


' ' 

u 

n 

= R 


u„ 


V 

s 

\ > 


, j 


To find the boundary values at the new time step we 
must calculate fourteen quantities corresponding to the 
seven dependent variables on each side of the interface. 
Since we are given four conditions we must supplement this 
by ten additional pieces of information. As discussed 
before we shall use the characteristic variables to ob- 
tain this information. We consider the interface as in- 
creasing in the counterclockwise direction. We denote 
by superscript 1 the region to the left and superscript 
2 the region to the right of the interface. The normal 
direction is taken as going from region one to region and 
the tangential direction counterclockwise. As stated 
before the signals c^, c^ propagating to the right bring 
two pieces of information to the boundary from the interior 
of region one. Similarly for region two we consider the 
signals ”^s travelling from right to left. In 

addition, we have three pieces of information on each side 
that travels along stream lines. Using the 
characteristic variables (see [33]) we find that 
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(6.10a) 


(6.10b) 


(6.10c) 


(6.10d) 


(6.10g) 


(6.10h) 


(6.10i) 


(6.10j) 
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(6.10k) 
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In these equations the bar on the right hand side indicates 
quantities that already have been computed. In all the 
examples presented they were computed using nearest 
point extrapolation from the interior of the appropriate 
domain (see also [40]). The subscripts n and s refer 
to the normal and tangential directions as given by (6.9). 
The sound speeds c^, c^, c were defined in (6.7), while 
Pq denotes a reference density usually taken as the density 
of the marker particle at the previous time step. Note, all 
material properties differ between the two sides of the 
interface. 

Solving the system of equations (6.10) we find that 


(6.11a] 
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with 


(6.11b) 




/p 


( 2 )^( 2 ) 

0 


Thus, the normal velocity is a weighted average of the 

extrapolated normal velocities from the two sides, the 

weight being the ratio of the acoustic impedances. In 

addition there is a correction term dependent on the 

difference in the extrapolated normal stresses. As 

usual T s= p - In a similar manner we have 

nn nn 


(6.11c) 
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Similarly we can solve for u_ and S__ on each side. 

S 5 S 

The complete recipe for dertermining the values assiqnrd to 
dependent variables on each side of the interface is given 
by the following algorithm. 
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(i) Extrapolate, from the interior to the boundary, 
the deviatoric stresses S, pressure p, internal 
energy e, and velocity components u, v on each 
side of the interface. 

(ii) Transform the stress and velocity components 
from the z ,r-coordinate system, to the n,s-coordi- 
nate system. 


(iiia) 

Use (6. 11c) to calculate t , then 

nn 


(iiib) 

Calculate pressure 

on each side using 




i = 1,2. 



^ nn nn ' 


= 0, 
ns 

(iv) 

Apply the slip free 

condition; = 

ns 


i = 1,2. 



(V) 

Transform the deviatoric stresses bade 

to the 


z ,r-coordinaLe system using (iv) and the extra- 
extrapolated values , i = 1,2. 

S 9 

(vi) Calculate u using (6.11a); then 

n 

(vii) Convert the velocity components to the z,r- 

coordinate system using (vi) and the 
extrapolated values i = 1,2. 

(viii) Calculate densities from, the equation of state 

(p ,e ^^^ ) using (iiib) and 
the extrapolated values i = 1,2. 

We note that the continuity of normal stress 
is not used to calculate the deviatorics but instead is used 
to compute p. This algorithm has been chosen so that the 
formulas are valid even if one or both materials are purely 
inviscid. In the limiting case of infinite density ratit 
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across the interface# i.e., ■*■ 0» the above 

algorithm coincides with the common driver-driven model 

and u » 

nn nn n n 

The method described above applies equally well to free 
surfaces if we interpret one of the domains to exhibit 
material properties such that all dependent variables are 
zero. For definiteness denote the vacuum side by super- 
script (2) ; then 

"nn - “ ' 

and the slip free condition is 

(6.10c') X „ • 0 . 

ns 

In many problems there exist regions with large 
gradients. These gradients need to be resolved in order to 
have an accurate solution. To prevent an excess of grids a 
stretching is introduced which is a function of only one 
coordinate transform. We denote the transformed variables 
by (a,S) . By use of the implicit function theorem it is 
possible to transform divergence free quations from the 
physical space (z,r) to the computational space (a,P) in 
such a manner that the new equations are still in diver- 
gence free form. A consequence of this is that shocks 
will be computed with the correct jump conditions in the 
computational plane- 

In most problems in ballistics the large gradients 
occur near the axis of symmetry. For these cases we used 
a transformation (see (160]). 
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(6.12a) 





with b » 1/(1 - d/a). This maps the interval (0,a) into 
(0,1) and d is the boundary layer thickness. For the 
acoustic problem discussed in the next section an alternate 
stretching was used. 


(6.12b) 


r = P 


1 - ae 


-b? 


This is inverted numerically to find $ ■ 8(r). For the 
acoustics problem the transform was used in both the z and 
r directions. In both cases one of the free parameters 
controls the size of the gradients while the other controls 
the size of the boundary layer effect induced by the trans- 
form. 

Another alternative is to choose new coordinate points, 

in some a priori manner. Cubic splines can then be used 
to construct the derivatives at The coordinate trans- 

formation generated in this manner need not be monotone. To 
enforce monotonicity one may need to alter slightly the 
so that no abrupt jumps occur. Even though this procedure is 
not automatic it allows the construction of flexible coordinate 
systems. Since the mesh generation is done only for each 
geometry the inconvenience is not great. The author has used 
this method to construct a mesh system for a generalization 
of (4.1) that is being used for optimization studies for the 
national transonic tunnel at NASA Langley Research Center. 
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Du« to the complicated geometry of the tunnel it uaa necessary 
to refine the mesh in several different sections of the domain; 
this made it difficult to use analytic transformations. The 
disadvantage of this method is more pronounced when one 
wishes to frequently change the mesh or change the geometry 
of the configuration. 

In other problems > as the shaped charge, transforma- 
tions are needed because of the shape of the domain. Even 
if a domain is rectangular if it is not properly aligned 
with the coordinate axes a large waste of grids will occur. 

In the shaped charge the initial shape of the liner is a 
rectangle rotated about 30** with respect to the axis of 
symmetry. The length of the liner is about 100 times 
larger than the width. For this case it was necessary to 
introduce the computational coordinates (a, 6) via a full 
two dimensional transform. Even though this complicates 
the equations it presents an enormous savings in both 
computer time and storage. Due to the moving boundary it 
is sometimes necessary to introduce time dependent trans- 
formations. In all these cases the coordinate transform is 
given analytically. Hence, it does not map the physical 
region into a perfect rectangle but rather into a computa- 
tional region which is in some sense reasonable. 
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The code described above, SMITE, has been applied to 
numerous problems in elasticity and plasticity. These 
include application to bars impacting on other bars or into 
plates as well as problems with shaped charges and explo- 
sives. All these applications involve highly distorted 
boundaries together with the interaction of many shocks. 

To illustrate this method and its range of applicabi- 
lity we present two examples. In the first case we consider 
an aluminum sphere impacting on a tungsten disc. The initial 
configuration is shown in figure 6.1a The configuration 
after 25 microseconds is shown in figure 6.1b. The sphere 
has flattened out and is extended normal to the axis of 
symmetry. The disc has also been indented. Obviously, 
large distortions have occurred in both materials. These 
graphs were obtaining the system (6.4) - (6.6) with the 
inclusion of plasticity. 

The first example illustrated the interaction of two 
metals. In the second example wc consider a gas-metal 
interface. The initial configuration is shown in figure 
6.2a. The larger region is a gas with ignition occuring 
at the origin. A blast wave propagates through the gas and 
impinges on the metal. The solution after 35 micro- 
seconds is shown in figure 6.2b. The gas has expanded 
outward while the metal has split into tv;o sections. 

The larger section is a slug which contains most of the 
mass. In addition, there is a jet region which moves 
rapidly to the right. Additional examples are given in 
{311 and [33). 
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7. RADIATION BOUNDARY CONDITIONS 

In many problems of interest one is interested in 
solving the equations in an infinite domain. For computa- 
tional expediency one needs to compute in a bounded domain. 
One possibility is to map the infinite domain onto a 
bounded one. However, in many circumstances this mapping 
can aggravate the situation especially if the solution is 
oscillatory at infinity or the mapping has a singularity 
(see [79]). An alternative possibility is to insert an 
artificial boundary and then impose boundary conditions 
on this surface to simulate an infinite domain, i.e. there 
should be no reflections from the boundary back into the 
domain. Unless certain restrictions are met this will in 
general not be possible 1851 . 

In general one cannot construct boundary conditions 
that give no reflections. Instead one wishes to consider 
conditions which are in some sense better. The notion of 
better can be defined in many ways. Some of them are 

(1) the reflections decrease rapidly as the position of 
the boundary goes to infinity 

(2) the reflections decrease for longer wave lengths 

(3) the reflections decrease as the incident wave 
approaches in a direction more normal to the boundary 

(4) the reflections are decreased so that the approach 
to a steady state is accelerated. 

One approach to ilccreasinq reflections is to intro- 
duce a viscosity near the boundary or to introduce a sponge 
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Iay«r [148). With this approach it is not clear 
lihat effect the boundary treatsusint has on the interior 
dynmics. in addition it is difficult to improve these 

if one wishes to further decrease the reflections. 
Conditions 2 and 3 were used by Engquist and Majda [53], 

[54] to construct an asymptotic set of nonreflecting 
conditions based on pseudodifferential operators. The 
higher order n^thods require Pade approximations for sta- 
bility. Rudy and Strikwerda [162] have constructed, by 
heuristic arguments, a radiation boundary condition based 
on (4) . 

Gustaffson and Kreiss [85] have shown that in general 
one can not construct nonreflecting boundary conditions 
unless the behavior of the solution is known in the neigh- 
borhood of infinity. We adopt their procedure and con- 
struct boundary conditions which are based on an asymptotic 
expansion of the solution valid for large distances. As 
with all asymptotic expansions we expect reasonable results 
even when the artificial boundary is brought in quite 
close to the region of interest. Extensive numerical tests 
indeed confirm that the domain of integration can be very 
constricted when one uses the higher order boundary operators. 

Specifically, we consider the linearized Euler equations 
in cylindrical coordinates. 


PVq + V 


(7.1) u^ + ^2 


'^t (vvq + p)j. = 


0,z ^3 


vvo,z - uv, 


n 


id»re (\Xq{z,x) , represents the mtM flow. 7he 

mean density is assismed ^constant and is scaled so that 
P© " ®Q “ 

In order for it to be feasible to Integrate (6.1} in 
a bounded domain it is necessary to assume that the n^an 
flow and the forcing terms decay as r or z go to infinity 
(see e.g. [85]). Hence for sufficiently large d (where 
^ *r *^2) (7.1) can be approximated by the wave syston 
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or since p is linearly related to p 


(7.3) - 7^p = 0 

2 

where 7 is the Laplacian in cylindrical coordinates. 

To find radiation conditions we first consider (7.3) 
in spherical coordinates and let d represent the spherical 
radius. It is Icnown [11] and [59] that p has a formal ex- 
pansion in terms of travelling waves 
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where e represents the angular dependence of p. 
One can then verify that 
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We note that the approxisiation to ^^*5) “ 0 is 

just what one would obtain from a one dimensional character- 
istic theory. More generally, let 
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It follows from (7.4) that 

<’-8> V = °(pfer. 

It is shown in [11] that boundary conditions B p « 0 

m 

all lead to well posed problems in the sense of Kreiss [81]. 
Furthermore, one can show that errors in the solution 
decrease as d goes to infinity. In particular for the 
first order approximation (7.6) we can show that the error 
p between the solutions in the bounded and infinite domains 
satisfies 


(7.9) 
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g^ dS 
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where g = 0(l/d ). In general for the mth order approxi- 
mation (7.8) we can show that 


^3 

tdier« g « oa/d®*^^) and I 1*1 1 » III •III* appropriate 
norms. 

The extensicm to cylindrical coordinates is straight- 
forward since 

(7.11) r = d cos 6 , d^ = r^ + 

2 

z ® d sin 6 , tan ® = ~ 

Hence, all derivatives with respect to d can be expressed 
as derivatives with respect to z and r. In fact one can 
show by induction that only derivatives tangential to the 
boundary need appear. 

The lowest order approximations in spherical coordi- 
nates are given by 
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For cylindrical symmetry we have by (7.11) that 

^^•^3) |§=|fcos9 = - [I? sinej 

Hence, we can replace (7.12) by 

Bj^p = (p - u sin 6 - V cos 3] + ^ = 0 


(7.14a) 
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In order to illustrate the advantage of even the 
first order boundary condition we present one example. 
(7.1) was Integrated using a mean flow obtained from 
experimental data. The source was taken as a single 
monopole along the axis of symmetry. In figure (7.1) 
we plot the pressure as a function of time for a fixed 
axial point, in figure (7.1a) the characteristic (or 
Sommerfeld) condition 


(7.15) 


l£ + lE = 

3t 3r 


0 


was used. The spurious reflections from the boundary 
are evident. In figure (7.1b) we present the solution using 
the first order boundary condition, Bj^p * 0. This 
solution no longer has any spurious reflections. These 
comparisons were made with the boundary at fixed distance. 
Varying the position of the artificial boundary it was 
found that, using the first order boundary condition, one 
could bring the boundary quite close to the sources with- 
out loss of accuracy. For distributed sources, moving 
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sourcos or a guadrapole it became neceaaary to use higher 
order boundary conditions. Further comparisons are pre- 
sented in UoJ and [ll]. 

7he problems in dynamic jet acoustics present severe 
difficulties for any numerical scheme. One wishes to find 
the acoustic pressure in the far field which necessitates 
many grid points. However, the grid must be constructed so 
that there is high resolution near the origin to resolve 
the sources and the mean flow. In addition one wishes to 
verify long term patterns of the pressure. This requires 
accuracy over many time steps. To have any chance of solv- 
ing this problem requires attention to the methods 
described in the past few sections. The use of high order 
it^thods is required to limit the number of mesh points. In 
addition the artificial boundaries must be brought in as 
far as possible to further limit the number of mesh points. 
This requires boundary conditions that severely restrict 
the reflection of waves from the artificial boundary. 
Careful attention to the stretching of the grid is also 
necessary. With all the above considerations a realistic 
two dimensional problem still requires about 35,000 mesh 
points for reasonable accuracy. This required the use of 
explicit methods that could be executed very efficiently 
on a pipeline computer. This problem shows that a success- 
ful code for a large scale problem requires an efficient 
algorithm and boundary conditions; a careful attention to 
the physics as well as to computer architecture. 
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For many applications one is interested in the fluid 

dyoasie equations linearised about a s»an state a uniform 

state which is nonsero at infinity. Let 

represent tiie uman state in the far field and define 
2 

c* » YP«/p*‘ We assum that is sero. The two 

diffiensional analog of (7.2) is 

“t + "-“x * i Px ■ “ 

00 

(7.16) + u^Vj^ + X p » 0 

^00 ^ 

Pt “ Px • 0 • 

This is equi^mlent to the convective wave equation 

(7.17) - o2(p^^+Pyy) - 0 . 


By a change of variables this can be transformed to the 
wave equation. We then transform the boundary conditions 
B^p - 0 to this coordinate system. When all is finished 
we transform back to the (x,y) system. The first order 
boundary condition for (7.16) then becomes, [13] 


(7.18) 
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2 2 2 
X + y and r 


“I 1 

c* - u* 


2 2 
X + y . 


Nhan 


this raduces to tie previous boundary condition. When u^ 
is larger than the outflow is supersonic and no boundary 

conditions can be specified at the artificial boundary. The 
boundary condition (7.18) has been used for the Navier-Stokes 
equations with a subsonic outflow. The use of (7.17) sub- 
stantially increased the rate of convergence to a steady 
state solution. The unknown linearised state (p^, u^, c^) 
was taken front the previous time step. A generalisation of 
(7.15) for nonlinear problems is considered by Hedstrom (94). 

It is also possible to generalise this theory to the 
three dimensional Helmholtz equation 


(7.18) 


A p + k p 


The analog of (7.4) is 


(7.19) 


-ikr 


* f^(e,0 


The radiation boundary conditions are then given by 
B_p - 0 with 


(7.20) 




m 


m 

n 

j»i 


f. 


ik + 


This procedure is ordered by first operating with j ■ 1 
and continuing until j » m. In [12] energy estimates are 
obtained for the error using (7.20). Numerical computations 


prttsonted in (12] denonstrate that for many caaaa on* n*ed 
only tak* ten meah pointa normal to the boundary, indepen- 
dent of k, in aolving (7.18). The boundary conditiona 
apply to the Laplace equation when k « 0 in (7.18) and 
(7.20). Extenaiona to tt^ diroensiona are alao conaidered 
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S. JU^PLICATICHIS 

In ttM previous ssetions ws have tried to stress the 
interplay bet«Men the nunerical schme and the physics of 
the ptobltm. In particular %ie strongly feel that there is 
no such thing as a universal or ultimate scheme. The 
variety of phen<»Bena described by time dependent partial 
differential equations demands a variety of n^anerical 
meth^s. In this section we shall descrlJde sc»m difficulties 
that arise in particular situations. Due to a shortage of 
space we can only sketch these difficulties with a brief 
description of possible remedies. References to the litera- 
ture will be given for a more extensive discussion of indi- 
vidual situations. 

A. Shocks 

When dealing with nonlinear systems the solution 
does not always remain snx>oth even when the initial data is 
SR^oth. Instead surfaces of discontinuities, called shocks 
arise. The solution is smooth on either side of the shock 
and across the shock the solution is governed by jump condi- 
tions. These shocks are an inherently nonlinear phenomena. 
Linear discontinuities such as contact discontinuities can 
also occur. For a survey of the analytical theory of shocks 
see (118). 

The standard convergence theorems for numerical 
schemes are based on the assumption that the solution is 
smooth (see (1561). When shocks occur it is well known that 
one can construct reasonable schemes that converge to a 


solution with ths wrong shock location (s.g. [116]# [209]). 

Otim sltsmstivs is to sns ths di££er«ntisl equations only 
in s»(x>th regions. The i^ock itsel£ is followed explicitly 
using the Rankine-Hugoraoit relations supplemented by <diar-> 
acteristic data (see [140]# [152], [166]). For coBg>lex 
flows with several intersecting shocks this is difficult. 
Furthermore, there is the additional difficulty of predicting 
the generation of shocks that do not exist initially. The 
introduction of viscous terms further con^licates the 
method. 

An alternative to fitting the shock is the so-called 
shock capturing method. Lax and Wendroff [116] showed that 
if the equations are written in divergence-free form 

(8.1) ^t ^ f = 0 

and if the numerical method also has this property then 
when the scheme converges it will give the correct shock 
speeds. Because of the ease of use, this method has dominated 
the computation of shocked flows. We note that the use of 
divergence-free form is not necessary. Instead one can integrate 
(8.1) as a quasi-linear system and compensate for this by 
the addition of terms that depend on the mesh (see [205i]). 

It is also well known that the use of the divergence free 
form is not sufficient to give the correct shock speed. 

Once one generalizes the definition of a solution to allow 
discontinuous solutions then the solution to (8.1) is no 
longer unique. One must demand additional constraints, e.g. 
that entropy increase across the shock, in order to have a 





82 


unique solution (see [118J for additional inforauition) . For 
first order systems necessary and sufficient conditions which 
guarantee entropy satisfying discrete shock profiles are 
given by Itajda and Ralston 1132]. 

The use of the difference scheme across the shock has 

several disadvantages. In some cases it has been found that 

non-physical rarefaction shocks can appear ([90], [125]) i.e., 

the entropy condition is violated. This generally happens 

at a sonic line or a stagnation point where an eigenvalue of 

the system becomes zero. All the known examples for model 

equations exhibit troubles only when the coefficient in the 

equation passes through zero. The usual second order schemes, 

Lax-Wendroff or MacCormack, have a dissipative mechanism 

3 f 

which vanishes when an eigenvalue of in (8.1) is zero. 

Hence, in these circumstances the scheme is effectively non- 
dissipative and so it is no surprise that difficulties can 
occur. Another difficulty of higher order method is that 
overshoots can occur in the neighborhood of shocks even when 
the shock speeds are calculated correctly. 

The standard cure to these difficulties is to add an 
artificial viscosity term to the scheme. This viscosity 
should be constructed so that it does not affect the accuracy 
of the scheme in the smooth portion of flow. It also should 
not vanish at the sonic lines or the stagnation points. In 
practical computations several viscosities have been suggested 
which work reasonably well [29], [49]. Majda and Osher 
[130], [133] have suggested some viscosities for which 
they can prove the convergence of the scheme to the correct 
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solution. Iiorat and Peyret [120] have shown that one can 
lessen the ioqpact of oscillations by choosing the 
correct variant of NacOoraack's method (3.5). This becomes 
difficult for complicatcid flow patterns. Lerat [121] 
discusses the addition of nonlinear correction terms to 
reduce the oscillations. These latter studies are based 
on the TOdified equation approach [100], [171], [194]. 

These corrections have been of a very specialized and 
problem dependent nature. 

It is known that monotone schemes have the property 
that overshoots do not occur and that they give the correct 
shock locations [90]. Unfortunately, linear monotone schemes 
are only first order accurate. Crandall and Majda [44] 
have considered generalizations of monotone schemes as well 
as extensions to several dimensions by using splitting methods. 
They demonstrate that splitting the equations into one 
dimensional portions can have strange effects on the 
shocked solution. The definition of monotonicity for 
systems of equations is not clear. 

An alternative to using monotone first order methods 
everywhere is to use these methods only near the shock 
and to use a higher order method in the smooth portion 
of the flow [88]. Rather than using different schemes 
in different regions it is easier to automatically com- 
bine these schemes using hybrid techniques [89]. The 
monotone methods produce excessive smoothing of the shock 
profile. Several nonlinear corrections have been suggested 
[18) » [91] » [191], [208] to prevent this smearing. 
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A silkier technique is to use the schen^ without 
the addition of any artificial viscosity. At the con^le- 
tioa of each step one Miters the solution to reiiMSve any 
oscillations. The sin^est way to accomplish this is with 
a Shuman filter [88] . l<et u represent the solution to 
any finite difference scheme at time t. We then define 
the corrected solution to be 

(8.2a) • 

with 


(8.2b) 


0 < < i . 


For (8.2) to be second order accurate, in space, we require 
that 9j = 0(Ax) in the smooth regions of the flow. (8.2) 
can be viewed as the second step of a splitting process which 
adds the viscous terms Ou ) . Hence, stability is ensured 
whenever the basic numerical algorithm is absolutely stable. 

A reasonable choice for 9 is [92] 


(8.2c) = a 


~it2~ °i-n~ *°1 


where 0 < a < i and o is a function of u. For the 
fluid dynamic equations one frequently chooses a as the 
density. 

The various techniques discussed all reduce oscilla- 
tions in the neighborhood of a shock. It is not clear 
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one iiMds to thaso oscillati<ms except 

for aesthetic reasons. 'Ssobs see»i to be problem dependent. 
F^r pn>blsaa ifith ca^wiMbmt it is isq^rati^^ to prevent 
oscillations nAich nay falsely trigger the conbustion 
process. For steady state calculations with separate 
shocks the oscillations probably do little harm to the 
total solution. For dynamic situations with interacting 
shocks the situation is less clear. 

Shocks are inherently stable and compressive phenomena. 
Hence, even without shock fitting the shock is smeared over 
only a few mesh points even for long periods of time. 
Hot^ever, contact discontinuities continue to spread in 
time. Hence, in muXtidimensions where coarse grids are 
necessary one cannot resolve contact discontinuities ever 
long periods of time. One can try to convert these dis- 
continuities into pseudo-shocks [91]. Alternatively one 
can use fitting techniques only for contact discontinuities. 
One such technique was described in section 6 for the inter- 
face between different materials. 

Numerical evidence indicates that if higher order 
methods are used in smooth parts of the flow then errors 
in the shock area do not propagate into the smooth region 
(188]. Since contact discontinuities are a linear phenomena 
it seems that errors in the discontinuity region may 
contaminate the entire domain of integration. For 
linear problems with discontinuities a straightforward 
method will yield only second order accuracy even in 


me^fSk xegiioiiMB (12H. 71^ aTOuracy of the aathod can be 
recovered by pre- and {^i^-processing {138}. For non- 
liiMar probloas preprocessing and postprocessing will not 
urorlc. On tile other hand there are indications that there 
is no need for any adjusti^nts for shocks (fll9], [188]). 

In fact preliminary c<»iputation8 demonstrate that one can 
achieve one point shocks and contact discontinuities 
using Chebyshev spectral methods (3.16 - 3.18). Small 
oscillations appear which can be removed by a postprocessor 
(0. Gottlieb, private communication) . 

When using implicit methods to compute the solution 
to shocked flows, computational experience indicates that 
one can not use time steps more than about three or four 
times the local Courant limit, even for stationary shocks 
(e.g. [93]). In many practical situations this is not a 
serious limitation as the time step for an explicit scheme 
would be goverened by regions other than the shock region. 
For example, in a shock -boundary layer interaction the time 
step is goverened by the boundary layer mesh and not the 
shock. Hence, an implicit method can still use time steps 
about fifty times larger than those used by an explicit 
method. 

B. Multidimensional Problems 

In several of the previous chapters we have dis- 
cussed methods for one dimensional problems. For practical 
applications one needs multidimensional codes. For some 
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nethodfl as leapfrog it is straightforward to instruct a 
stable oultidinensional ^benie. Even in this case the most 
straightforward version can lead to restrictive time steps. 
Consider the equation 

(8.3) BUy = 0 . 

The straightforward leapfrog method for this equation 
is 


(8.4) 
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This is stable if 

(8.5a) A sin S + B sin n _< 1 

for all C and n. In the worst case this can demand 
1^11 1 if I! Bji £ i. For the fluid dynamic equations (8.4) 


is stable when 


(8.5b) 


At 


1 fjLl 

Y(axJ iAyJ 



< 1 . 


Abarbanel and Gottlieb I 2 ] have pointed out that one can 
improve this stable condition by averaging the derivatives. 
In particular we can replace (8.4) by 


it 


( 8 . 6 ) 


,n+l 

‘i. 





^ on 

2Ay ®i,j 

j-l" 

This is now stable when 

(8.7) A 1 ^ » "iyS 1 1 

which is optimal. Similar extension to three dimensions 
exist except that the averaged scheme is no longer optimal 
but still better than the standard leapfrog method. Exten- 
sions to arbitrary dimensions are considered in [201]. 

A similar gain is achieved for the fourth order leap- 
frog method (S. Abarbanel and D. Gottlieb, private commu- 
nication) . 

For implicit methods the extension to multidimensions 
is less straightforward. The obvious generalization of one 
dimensional schemes leads to the necessity of inverting 
large sparse matrices which are no longer tridiagonal. 

To avoid this problem one usually employs an alternating 
direction method to reduce the problem, to a sequence of one 
dimensional problems. WcDonald and Briley (128] have stressed 
the importance of doing this in such a mianner that each 
portion of the split is consistent with the original equations. 
In chapter 5 we have indicated some of the disadvantages of 
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A.O.I. Thestt include a reduction in the rate of convergence 
to a steady state and also instability for some important 
three dimensional versia»as. 

For the usual multistep explicit methods one can con- 
struct full two dimensional versions. This has been done 
by Richtmyer [156] ,Burstein [29], MacCormack [124], Turkel, 
Abarbanel and Gottlieb [186] among others. These schemes 
usually have a restricted stability criterion. More impor- 
tant, it is difficult to treat different directions in a 
different manner in order to take advantage of the physics 
of the situation. Strang [179], Yanenko [206], and 
Marchuk [134] have introduced the concept of splitting 
the eqxiation into several components. We consider the 
general equation 

(8.8) ^ 4 . ^ M ^ M u ■ Mu . 

w X y z 

Here, we have arbitrarily split the right hand side into 
three portions. Frequently these splits are identified 
with separate dimensions. However, in other applications 
other splits are indicated by the physics (see e.g. [135] 
and section D) . We now consider the subsystems 


(8.9) 


V. = M v 
t X 


MyV 


t V, 


M V 
z 


and we denote the numerical solution to these subsystems by 


(8.10) * L v” ; = L v" ; = L v” 

X y 2 
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rcm^cti^ly. We can reconstruct the solution to f8.8) by 
(8.11a) u”'*’^ - L L L u” . 

X ^ z 

This is a first order, in time, approximation to (8.8). 

In order to make the approximation second order in time we 
follow (8.11a) by 

(8.11b) - L L . 

X y z 

This is stable if each of the one dimensional operators 
are strongly stable. Gourlay and Morris [75], [76] have 
considered the implementation using multistep methods. 

Gottlieb [67] has shown that this 2-cycle of permutations 
is second order even for nonlinear equations. This splitting 
in two or three space dimensions coupled with some multi- 
step scheme has been very »accessful for many applications. 

One can show that one achieves the optimal time steps. One 
can also treat different directions in different manners. 

For example, MacCoriaack [125] uses the operator in the y 
direction more of tv a, with smaller time steps, than the 
operator in the x direction. This compensates for the 
finer mesh in the y direction. Alternatively, one could 
couple Fourier methods in periodic directions with finite 
difference or Chebyshev methods in directions with boundaries. 
For a further description of details see [77]. 
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As with all awttods splitting techniguss hava their 
drai^cka. One inaediate difficulty is that the intenimdiate 
steps have no physical interpretations. Hence# if the 
coefficients of the equations or the boundary conditions 
depend explicitly on time it is not clear which time to 
use for the middle steps. This is compounded if multi- 
step methods are used for the one-dimensional operator. 
Additional difficulties are encountered wh*?n the solution 
has a shock inclined to one of the coordinate directions. 
Crandall and Majda [45] have shown that unusual occurences 
happen in this case. The efficiency of the random choice 
method also deteriorates in the presence of oblique shocks 
when splitting is used (Chorin, personal communication) . 

Another difficulty occurs when some subsidiary con- 
straints are intrinsic to the solution. For example# for the 
Maxwell equations the condition div B » 0 can be viewed 
as an initial condition. For problems with variable 
coefficients (8.10) is not symmetric in x and y even 
though the operators appear in a symmetric fashion. Hence# 
numerically div B iu nonzero. This introduces a non- 
physical Lorenz force which can cause numerical instabilities. 
The author has done extensive calculations with the nonlinear 
ideal MHD equations, ifhen dimensional splitting was used 
nonphysical instabilities always occurred. This also 
occurs when Lagrangian schemes are used in MHD (Brackbill 
and Barnes [21]). A similar situation occurs in the 
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in^»npresiibl« Navler-Stoktts equations tritere divu • 0 
nust be satisfied. In this case the situation is inproved 
by explicitly adding tezM involving divu to the 
equations. This does m>t affect the analytic solution 
since divu - 0. However, the numrical solution can be 
stabilized by such a procedure (see (86], (99]}. 

C. Aerodynamics 

Aerodynamics is frequently divided into internal 
and external ilowB. Internal aerodynamics describes flows 
in nozzles, ducts and turbomachinery. These flows usually 
involve complex flow patterns and frequently require the 
addition of chemistry models to study the propagation of 
flames. Boundary conditions are extremely important for 
internal flows. These flows are usually subsonic or 
transonic with maximum Mach numbers of about 1.3. Hence, 
one is frequently not Interested in shocks and can frequently 
dispense with the conservation forms of the difference 
equations. This is especially important for implicit methods. 
The block structure of the matrices to be inverted are much 
more complicated when the divergence form of the equations 
must be used. The velocity form of the equations frequently 
allows the decoupling of the blocks into a direct sum of 
smaller blocks (see (23]). It is much faster to invert 
three scalar tridiagonal matrices than to invert one block 
tridiagonal matrix with 3x3 blocks. Hence, it is computa- 
tionally important to utilize the proper form of the equations. 
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Vor myutmm fr<Mi chmicAl notela this !■ svsn 

MOTS i^)ortsiit. Ah 6xt«Hii^ discussion of ii^licit mthods 
for intsrnsl flows is giswii by McDonald and Brllay 1128]. 

Por both intamal and sxtarnal flm^s ona must usa a nosh 
tiiat rssoli^s the boundary layero. This can be done by 
a nesh stretching or a finite voIism technique. In nany 
situations one is only interested in the steady atate solution. 
For this case iic^licit nethods are becking popular. As 
nsntioned in section 5 the i^>licit nethods becone nore 
inefficient as the complexity of the equations increase. 

This complexity can be created the existence of nany 
equations or the existence of cross derivatives. For nany 
aerodynanic flows one simplifies the full Navier-Stokes 
equations to the so called "thin layer" approximations [175] 
so that the flow is essentially unidirectional [51]. This 
is useful only when the boundary layer has some simple 
structure parallel to the body. An alternative is the 
rapid solver proposed by MacCormack [126]. Details of these 
approaches are given in the article by Hollanders and 
Vi viand [101]. An alternative to using the time dependent 
equations is trying to solve the steady state equations 
directly. Iterative methods have not proved to date, to 
be very promising. One possibility is to use Newton's 
method for linearizing the problem coupled with Gaussian 
eliminatio' . This has been carried out by Blomster and 
Skdllenw) [17] and Rizzi []57]. The major difficulty 
with this approach is that the bandwidth of the matrix 
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inoTMMS rapidly aa finar iMshas ara naadad for mora 
i^liatic aituationa. ^lia ia aafwcially aavara in thraa 
apaca diawnaiona. Altaraata linariaationa, auch aa 
8chubart*a natlK>da# aiay allaviata tha difficulty ainca ona 
can account for tim aparaity of the matrix to ba inverted 
(aaa a. 9 . [IHU (167]}. 

An iB^rtant problem in uaing the time dependant 
aquationa to achieve a ataady atata it finding waya to 
accalarata tha convergence to the ataady atata. One auch 
matl^d ia to uaa different time atepa at different R»ah 
pointa. For explicit mathoda tha time atepa could ba 
choaan ao aa to aatiafy tha local atability limit (Buratein, 
paraonal communication} . For implicit nwthods large tiiM 
atepa can ba viewed aa an iteration parameter. Some ideaa 
for optimizing thia parameter were diacuaaed in chapter 5. 

Even more important ia the proper implementation of 
initial and boundary conditiona. Rudy and Strikwerda [163] 
have conaiderad the effect of varioua boundary treatmenta 
on the acceleration to ateady atate. Over-apecificaticn at 
inflow can accelerate the convergence but frequently leada 
to oacillationa in the ateady atate. Under 'apecification 
at outflow can prevent the achievement of the ateady atate. 
The uae of a radiation boundary condition can dramatically 
decrease the number of atepa required to achieve the 
ateady state. An analytic treatment of several boundary 
treatments is provided in [147]. 
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The cofrect choice of ini^al 4ete ip elso important 
in achieving a rapid steady state* f^en the initial data 
does not satisfy the boundary conditions waves can arise 
idiich take a long time to dissipate. One procedure is 
to use a primiti\^ form of the multi-grid TOthod [22] . 

In this case ve solve the equations on a coarse grid, 
nils can be d<»ie rapidly and so the initial guess is not 
iiqiortant. The converged solution is then used as an 
initial guess for the next finer grid. A pyramid of 
grids can be used until the finest grid is reached. Most 
of the work is done on the finest grid and so the 
coarser girds do not increase the work by very much. By 
providing an excellent initial guess this process can sub- 
stantially accelerate the process of reaching a steady 
state . 

For complex three dimensional flows it is likely that 
all the processes described above and in chapter 5 will 
be necessary to achieve a steady state within a reasonable 
number of iterations. For surveys of methods for the 
Navier-Stokes equations see [7], [35], [137], [154]. 

In addition there will be, in the future, more of a 
need for internal accuracy checks, e.g. changing the 
mesh size. This is crucial for investigating the effect 
of various acceleration techniques on the accuracy of 
the steady state solution. In most studies the only 
accuracy checks are comparisons of averaged quantitities 
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With •xpwrimenfal aata, #.f. skin friction, drag, etc. 

These comparisons ignore any oscillations or other diffi- 
coll;tes that toy occur locally in sections of the flow. 

d. Meterorology 

The problems that occur in computational meteorology 
are very different from those in aerodynamics. First, 
in toteorology one is interested in the dynamic behavior 
of the model. Even in climatology one wishes to find 
statistical averages rather than a steady state solution. 
Furthermore, the equations of motion are a classical 
example of a system with different time scales. Allowable 
speeds in the atmosphere range from sound waves down to 
Rossby waves. The physically dominating Rossby waves 
are about 20-30 times slower than the speed of sound. 
Nevertheless, if one uses an explicit scheme the time 
steps is restricted by the fastest possible modes of 
propagation. A further restriction on numerical methods 
is that one wishes to limit the amount of dissipation 
introduced by the scheme. The long term weather patterns 
are governed by a delicate balance between heating due 
to solar radiation and dissipation due to friction and 
interactions with the oceans. For long term weather pre- 
diction it is essential that the numerical method does 
not interfere with this balance. 

As a simplified set of equations we analyze the 


shallow water equations in cartesian coordinates. 


( 8 . 1241 ) , 4 - VUy * 

(8.12b) + uVjj + Wy - -ghy - fu 

(8.12c) + uhjj + vhy = -h(Ujj + Vy). 


(u#v) are the velocity components and h is the equi- 
valent height of the atmosphere, g and f are ta)cen 
to be constants. The phase speeds of the system are 


(8.13) 


= u sin e + V cos 9 



In this simple system o)^ represents the important Rossby 
wave while <^2 3 relatively less important 

gravity waves. The flow is called geostrophic if the 
right hand side of (8.12a) and (8.12b) are zero. The 
flow is incompressible if the right hand side of (8.12c) 
is zero. The real atmosphere is quasi-geostrophic and 
aln^st incompressible. This is the cause of the small 
amount of energy in the gravity waves. 

One way to ovej.’come the difficulty with the different 
time scales is to use a semi-implicit method. In this 
method the right hand side of (8.12) is treated implicitly 
while the convective terms are treated explicitly. The 
resulting stability condition then depends only on the 
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velocities and not on the sound speed, A farther 

analysis of this is presented in Elvius and Sundstrom [52). 
Kavon [144] has ccnabined a semi-implicit method together 
with the linearization technique described in section 5. 

He found that it was necessary to iterate the procedure 
once in order to maintain the accuracy. Isaacson, 

Marchesin and Zwas [104] have used a fully implicit compact 
fourth order method coupled with linearization algorithm. 

The 8chen« (3.9) was generalized to two dimensions by using 
an alternating direction method. They found that the 
correct treatment of the singularity at the pole was 
necessary to maintain stability. The order of the factors 
in the A.D.I. method was also crucial for stability. 
Williamson [200] has compared the effectiveness of high 
order schemes for the primitive equations. He found that 
the horizontal diffusion term had a greater effect than 
the order of the scheme. The introduction of topography 
further complicates the comparison since for realistic 
grids large mountain chains are described by relativity 
few grid points leading to large gradients in the vertical 
variable. 

Gadd [62] has suggested splitting the system into fast 
and slow components and treating each separately with a 
Lax-Kendroff method. This has the disadvantage of using a 
dissipative method. It is possible to use a similar idea 
with the leapfrog method by using two grids. The convective 


temt are dlfferezuired on the fine grid tdiile the right 
hand side of (8.12) is differenced on a coarse grid but 
with a contact fourth order method. The resulting 
approximation to (8.12) is given by 
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An additional difficulty that occurs in tneterology 
is due to the spherical coordinate system used for the 
globe. The coordinate singularity at the pole together 
with the convergence of the latitude lines at the poles 
forces an unrealistically small time step for explicit 
schemes. Many attempts at using other coordinate systems 
or patching several coordinate systems have not been very 
successful. To enable the use of larger time steps some 
smoothing algorithm is used near the poles to eliminate 
the higher frequencies. In figure (6.1) we display a 
contour map of h for the solution to the shallow water 
equations in spherical coordinates. This solution used 
3 minute time steps and includes Fourier filtering near 
the poles (see [190], [199]). In figure (6.2) is shown 
a second graph using the scheme (8.14) together with 
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Ponri«r filt«rin9* This allowed 7| mintite tine steps. 

A leap frog method without any filtering could use only 
I minute time steps with the grid used. 

More realistic models are based on the primitive 
equations and also contained many levels in the vertical 
direction. Usually a pressure-like coordinate is used 
in the vertical direction. Even though many calculations 
have been performed with the primitive equations 
Oliger and Sundstrom [147] have shown that the system is 
not well-posed. Browning [26] has proposed a substitute 
which is a proper limit of hyperbolic equations. 

In volume 17 of Mztkodi oi Computational Pktjiia a 
variety of methods are presented for solving these equations. 
Spectral, pseudospectral and finite element methods all 
present advantages and disadvantages for large scale 
problems. To date there has been little comparison of 
these numerical techniques under real-life situations. 

e. Combustion 

The replacement of an ideal gas by a real gas 
and the inclusion of chemical processes introduces 
several new difficulties into the computation. Because 
of the exponential dependence of the ignition on the 
termperature it is critical that there be no overshoots 
in the computation of shocks. False overshoots can ig- 
nite the chemistry at entirely incorrect places and times 
and completely invalidate the computation. As previously 
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mentioned one alternatl^ is to introduce severe dissipa- 
tion. Anti-diffusion is later added to sharpen the 
shock profile (see [18] > [91]# [191]). An additional 
difficulty is the Introduction of many new equations 
when many chemical species are present. This is a parti- 
cular difficulty for implicit methods as the work increases 
with m^ where m is the number of equations. In certain 
cases Briley and McDonald [24] have shown that these 
matrices can be simplified. 

Since the flame occupies only a small region of the 
computational domain it is necessary to introduce a fine 
mesh which moves adaptively in time. Efficient ways of 
changing the grid, especially for several space dimensions 
are unknown. One way of moving the mesh for one space 
dimension is given in [50] • An additional difficulty with 
combustion problems is the stiffness of the ordinary 
differential equations which describe the chemistry. This 
necessitates much smaller time steps for the chemistry than 
for the fluid dynamics. This can be partially alleviated 
by using splitting techniques (section b) to split the 
chemical and fluid dynamical portions of the calculation. 

A different approach to these problems was introduced 
by Chorin [38]. He extended a probabilistic method due 

to Glimm [65] into a practical method. The random choice 
method is based on ,ne solution to many Riemann problems. 

In each interval the solution is considered to consist of 
two constant states with a discontinuity between them. The 
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solntion is advanced to tha next time level by analytically 
solving this local Rienaan problem and then choosing the 
new constant state to be the solution evaluated at a randcxn- 
ly chosen point within the mesh. This ^thod has the pro- 
perty that the discontinuity is sharp (in a probabilistic 
sense) without overshoots. Hence, there is no need for 
artificial dissipation or a moving grid. A drawback of the 
method is that it is at best first order accurate in the 
smooth regions of the flow. Since convergence is guaranteed 
only in a probabilistic sense accuracy for any given compu- 
tation may be poor. This depends crucially on the choice 
of sampling. It is expensive to calculate exactly the solu- 
tion to a Riemann problem at every mesh point and each time 
level. An alternative is to solve the Riemann problem with 
a finite difference method [93]. This also has the advantage 
that it easily generalizes to complex systems of equations. 

The present extension to several dimensions it. based on 
splitting techniques. As this reduces the effectiveness of 
the method other approaches are being investigated [66]. 

Sod [173] has compared this method with some finite difference 
schemes for simple shock problems. 
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f. PlaiRia Physics 

Ths simplest set of equations that describe a 
plasma are the ideal MHD equations developed by Lunquist 
1123]. These form a set of eight nonlinear symmetric hyper- 
bolic equations. Grad [78] has shovm that the full equations 
include all of standard fluid dynamics within a small range 
of parameter space in the MHD equations. Hence# all the 
difficulties discussed in the previous sections automatically 
occur in MHD. In order to make the problems manageable 
one can only treat simpler problems than those solved 
for the Navier-Stokes equations. This is even more true 
when more relevant equations# as the Vlasov or Fokker- 
Planck equations are considered. For these equations 
only one dimensional or simple geometries can be consi- 
dered. 

As a simple illustration we consider the steady 
state equations. For the Navier-Stokes equations we can 
solve the full system with the time derivatives set equal 
to zero. Some of these methods were described in section 
c and others are described in more detail by Hollanders and 
Viviand [101]. For ideal MHD even the case with no flow# 
i.e.# u ■ 0# is nontrivial. In this case the equations 
become 


7p =7><7xB 


(8.15) 


div B 
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syst«i of four oquations for p, By, B^ hat two 
real and two iaaginary characterittics. Hence, this it 
a adxed hyperbolic-elliptic tytten. For three dimentional 
aymtwm in a torua it it not even clear that the differential 
eyaten it well poted. By impoaing additional conatrainta 
Betancourt and Garabedian [16], [9] have* obtained aolutiona to 
(8. IS) for three dimentional configurationa. Brackbill {20] 
haa obtained aolutiona by uaing the full time dependent MHO 
equationa and marching to a ateady atate while removing 
kinetic energy to that u ■ 0 in the limit. It it 
o' /loualy much more difficult to do the type of ateady 
atate calculationa that are c^n^nplace in aerodynamic 
flowa. 

In addition to the ateady atate, or equilibrium 

problema, one ia alao interested in the atability of the 

equilibrium. When one wishes to study nonlinear stability the 

main technique is to integrate the time dependent equations. 

As n^ntioned in section b one difficulty is introduced by 

a 

the constraint div B *0. If divB is nonzero 
numerically it can introduce false sources of instability. 
Hence, stable equilibria can appear as unstable equilibria 
(see also {21], {142]}. 

For the time dependent equations there are three sound 
speeds known as the fast (or magneto-acoustic) speed, the 
Alfven speed and the slow speed (see e.g. {39]). Many of the 
instabilities in a torus occur at the Alfven speed which can be 
much smaller than the fast speed. Hence, we again face the 
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{^•rnnwiia of difforont tiiM ■calos. W« wish to use tia» 

■tops basod on the Alfven speed without sacrificing stability. 
One approach is to use an implicit, alternating direction 
method (see {1221). However, this requires the inversion of 
block 8x8 tridiagonal matrices which is tiiM consuming. 

Plasma flow has the property that s»st of the change occurs 
along magnetic flux surfaces and not across them. By using 
flux surfaces on a coordinate system one can separate out 
the fast and slower motions. Since the Alfven wave is 
incompressible this can be used to set up a scheme which 
uses tiiM steps based on the Alfven speed (see (106]}. The 
main disadvantage is that the method fails when the flux 
lines no longer provide a reasonable coordinate system. 

Another approach is to use a semi-implicit method (Brackbill, 
priva;.e communication) similar to that in section d. One 
can now use a standard Eulerian or Lagrangian grid. The 
si^e of the matrices to be inverted are severely reduced 
while the tinw step is governed by the Alfven time step. 

When considering non-ideal effects the major effect 
is resistivity rather than viscosity. While this still 
Introduces parabolic terms nevertheless the physical 
effects are quite different. It is well known that 
resistivity can frequently cause instabilities that did 
not exist in the ideal case [48]. As with the Navier-Stokes 
it is necessary to resoi w boundary layers where the 
resistive effects are dominant. An additional difficul ty 
with plasma physics is the existence of a vacuum outside 
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I th« planu. Th« aagnctic field in the vecuun ie governed 

by an elliptic equation which nust be couped with the time 
dependent MHD equationa across a moving interface. When 
solving the plasna-vacuvsn problem Lagrangian type methods 
or flux surface methods have been used. The Eulerian 
Mthods usuiliy represent the plasma as a low density plasma. 

We have concentrated on the application of MHD to 
magnetic controlled thermonuclear reactions. Another 
application is to astrophysics. For these problems 
radiation boundary conditions as described in chapter 
7 are i^ortant. Another important area is laser fusion. 

In this area the introduction of artificial boundaries 
is also in^rtant. In addition, these problems frequently 
have regions with extremely different properties as the 
pellet is compressed. Hence, the methods must be capable 
F of handling changes in coefficients in the range of 15 

orders of magnitude. Further discussion of methods for 
plasma physics are given in (155] and also volume 16 of 

f 

kiethodi oi Computational PlujsicA. 

I 

I 

f 

t 

g. Other Applications 

In this brief survey we have shown some of the 
difficulties that occur in specific applications. 

Naturally this survey cannot cover all topics. One major 
field which has not been discussed is two phase flow. 

Aspects of that field are discussed by H. Wirz [203]. 



Another field %ihich was touched on briefly in chapter 8 
is acoustics. This includes such diverse fields as 
seismology, underwater acoustics and jet acoustics. Each 
of these present their own difficulties. In seismology 
in particular one tends to use second order equations 
rather than first order systems. The occurrence of 
layered media introduces other difficulties. Both interior 
atr^ exterior ballistics provide problems with extreme dis- 
tortion of materials. Due to the high pressures and tem- 
peratures even metals deform. Hence, one cannot 
integrate in a region with fixed boundaries. Instead, 
the motion of the free surfaces and interfaces must be 
calculated. An additional difficulty is the existence 
of a plastic regime. For this regime the deviatoric 
stresses remain on the yield surface, i.e., 

j = K. It is a nontrivial problem to enforce this 
constraint numerically. A survey of numerical methods 
in elastodynamics is given in [5]. 

One major field which has not been discussed is that 
of incompressible flow. Discussion of finite difference 
methods is given in a companion article by Krause [109] 

(see also [61] and [185]). There has also been much work 
in the application of finite element approaches to solving 
the time dependent incompressible Navier-Stokes equations 
(see e.g. [58], [25], [102]. Of particular note is the 
recent use of grid-free methods. Chorin [37] used a method 
based on the interaction of a finite number of vortices. 
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Coabination grid-vortex methods have considered in 1271, 
U27], and C167J. Peskin [1531 has advocated the use of 
a Lagrangian grid free method. These codes have been used 
for many aeronautical purposes. The use of these codes for 
biological studies is complicated by the fact the boundaries 
which represent tissue material, are permeable [1511. The 
spectral methods can also be viewed as grid-free methods. 
Applications of spectral methods to incompressible flows 
are surveyed in [150] 

Boundary layer computations have been dominated by 
the use of fourth order finite difference and finite 
element algorithms (e.g. [41], [107], [161], [196], [205]). 

A major difficulty in the practical use of these codes is 
the lack of sufficiently accurate turbulence models. 
Extensions of the mathematical models to detached and 
reverse flow is also being investigated [34]. The boundary 
layer equations are a parabolic system. Other parabolic 
problems include heat flow and diffusion problems. Diffu- 
sion problems are of major importance fcr such diverse 
fields as oil studies, plasma diffusion and biological 
processes. 

Transonic flow calculations have been mainly based on 
the potential equation. Recent progress in the field has 
been presented in [204] and in the Proceedings of the 
Fourth AIAA Computational Fluid Dynamics Conference. 
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